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This  collection  includes  work  with  numerical 
hydrodynamic  methods  of  weather  forecasting.  A 
number  of  articles  is  dedicated  to  certain 
aspects  of  objective  analysis  of  meteorological 
fields  and  to  development  of  methodology  of 
operational  forecasting  by  numerical  methods  with 
the  aid  of  electronic  computers. 

The  collection  is  dedicated  to  scientific 
colleagues,  engineers,  weather  forecasters, 
students  of  higher  courses  in  hydrome  eorologic 
Institutes  and  universities,  interested  in  the 
problems  of  calculating  weather  ahead  of  time  by 
hydrodynamic  methods. 
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IMPROVEMENT  OP  THE  METHODOLOGY  OP  FORECASTING 
THE  BARIC  FIELD  FOR  SEVERAL  DAYS 

S.  A.  Mashkovich 


A  generalization  is  given  of  an  earlier  proposed 
[2]  linear  three-level  forecast  model.  The  effect 
of  horizontal  turbulent  mixing  in  the  vorticity 
equation  and  effects  of  surface  friction  are 
considered  additionally.  The  system  is  com¬ 
puterized.  An  evaluation  of  the  forecasts  is 
given. 


Works  [2],  [3]  proposed  a  forecasting  system,  which  is  based 
on  a  linear  three-level  model  of  the  atmosphere.  The  formulation 
of  the  problem  was  based  on  a  method  proposed  by  Ye.  N.  Blinova. 

This  system  is  used  in  the  work  of  the  hydrometeorologic  center  of 
the  USSR  for  forecasting  ground  pressure  for  several  days. 

The  present  article  is  dedicated  to  the  improvement  of  the 
mentioned  forecasting  methodology.  Namely,  in  the  vorticity  equation 
a  component  describing  horizontal  turbulent  perturbation  is  con¬ 
sidered  additionally;  furthermore,  the  effect  of  surface  friction 
is  drawn  in. 


In  this  Instance  linearized  vorticity  equations  and  heat 
inflow  equations  are  written  respectively  in  the  form: 


FTD-MT-2^-174-70  1 


2»«3fco«l  dpm  ,  k. 

—p - 3T+7aa»- 


\*t  »x)  a:  rft  a>  2m cot i  c 


-J- 

■5 


Ail. 

a: 


(1) 


!!ere  the  following  designations  are  introduced:  -  stream  :function, 

u  -  vertical  component  of  velocity,  t  -  time,  0,  X  -  spherical 
coordinates,  u>  -  angular  velocity  of  earth  rotation,  a  -  index  of 
circulation,  C  ■  p/P,  p  -  pressure,  P  -  1000  mbar,  aQ  -  radius  of 
earth,  g  -  acceleration  of  gravity,  and  feg  -  coeff.  cients  of 
turbulent  mixing,  R  -  gas  constant,  y  and  y  -  dry  adiabatic  and 
vertical  temperature  gradients,  A  -  Laplacean  operator  in  spherical 
coordinates . 


Let  us  introduce  dimensionless  variables  according  to  the 
formulas 


t  — tt\  p-pV. 

.t-  - 


and  designate 


i-S-i)  Pft  i  -  * 

- - - ~ - T~-  *■“ 


Then  system  of  equations  (1)  takes  the  form  (primes  on  letters 
have  been  omitted): 


-£i_  + ; -21±_ 
ditfc  a><>: 


jLl  A5: 

4*  0'r> 


Tip* 


1  dt  4) 1  ■  m  AT 


A 

4- 


4ii 

AC  ‘ 
■AA^. 


(2) 


We  introduce  the  following  designations  for  the  operators: 


z-  £+’aT-  it1' 


(3) 
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Then  equation  (2)  can  be  written  in  the  form 


.0 

JT 


(4a) 
(4b ) 


We  consider  processes  in  the  lower  half  of  the  atmosphere, 
namely:  in  the  lajer  located  between  sea  level  (c  -  1)  and  the 
isobaric  surface  of  500  mbar  (c  «  0.5). 


As  boundary  conditions  we  take: 


when  t— 1 

(5a) 

—  d  when  i  — 0.5. 

(5b) 

Condition  (5a)  speaks  about  the  fact  that  at  the  lower  boundary  of 

the  considered  layer  of  atmosphere  exist  vertical  motions  conditioned 

by  the  effect  of  friction  in  the  boundary  layer  (see  [1]).  Condition 

(5b)  designates  assumption  about  the  quasi-barotropic  nature  of 

motions  in  the  middle  part  of  the  atmosphere.  It  allows  writing 

equation  (4b)  at  a  level  c  ■  c  in  the  form 

c 

/.•i  —  0  when  :  — (6) 


Consequently,  the  stream  function  at  a  level  C  ■  C  can  be 

c 

found  from  equation  (6)  independently  of  the  solution  or  system  (4). 
The  solution  of  the  equations  (6)  is  written  in  the  form 


where 


V* 


Re  '  *•"  ’/?(*). 


v»  +  (yf;  *• 


y  »  II. 


(7) 


P™(9)  -  associated  Legendre  polynomial. 
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Let  us  break  the  considered  area  of  the  atmospheie  into  two 
layers  (see  Fig.  1).  We  will  consider  a  three-level  model  of  the 
atmosphere  including  sea  level  (s  *  l),  500  mbar  (c  *  t  )  and  a 

w 

certain  intermediate  level  c  =  ;  . 

n 


Fig.  1. 


Integrate  each  of  equations  (4)  over  the  vertical  coordinate 

within  limits:  1)  from  ;  =  1  to  ?  =  C  ,  2)  from  C  =  t,  to  ?  a  t  . 

n  n  o 

Here  indices  1,  n,  c,  mean  that  the  quantity  is  referred  respectively 
to  sea  level,  intermediate  and  mean  levels.  When  the  integrals  are 
not  calculated  analytically,  we  will  approximate  them  by  method  of 
trapeziums.  Then  system  (iJ)  with  allowance  for  boundary  conditions 
(5)  can  be  written  in  the  form: 


where 


7  •  7  •  r>  T  ^  '^i  in 

*15*1  “  ^-nVn  —  *A,at  2 


rnwn  -f  r,*A 

: - 2 — — «,-=o. 


/  7  \  -  9A  n  **  '**'  I'nWn  -t-  r,«-r  n 

*•*■*  ox  2 - f-77«»-0. 

+  C'in)  f-i^i  ’>i  —  g*«*B  *=  0. 


—  P<«V  “r  0, 

“  jji*'  ~  ?«)•  ~  *)• 

fll  **  I  —  <*2  **=  —  S*  . 


(8) 


System  of  equations  (8)  is  heterogeneous  and  contains  a  known 
right  side  of  type  (7). 


Let  us  find  the  generax  solution  of  uniform  equations,  which 
correspond  to  system  (8)  We  look  for  this  solution  in  the  form: 


m,  m 

l(ml  +  9*  I  \ 

*' “**«?/“■*  ]pc(H)‘  (9) 

(7—1,  n,  c). 

To  reduce  the  writing  let  us  drop  subsequently  the  indices 
n  and  m.  We  introduce  the  designations: 

yn), 

z>  -  iy{9  -f  /  lyin)  -f  2/m  -  y*j. 

Substitute  (9)  in  the  system  of  uniform  equations.  After 
simple  transformations  we  go  to  the  following  equations: 

B‘P>  ~  a, A, ml  +  % r,*y]  _  b\{ZB  +  a&ml)  -  «  0; 

Bm[Za  —  a^Atmi)  —  Am  r,  —  Ac-~*~  *■*  b*{z,.-\- a^A^nl ) ; 

Bn-£-Ln  P^C  — 0,  (10) 

Equate  to  zero  the  determinant  of  system  (10) 

j  Z,  —  1C,  -f 
!  n 
>  C^Lx-\-Cu 
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Condition  (11)  gives  us  the  characteristic  equation  for 
determination  of  a,  which  can  be  written  in  the  form 

+  fc)o  +  </  +  //-0.  (12) 

where 

u  “  —  kiy*  -f  +  3*.i)y  —  Pm; 

b  =  -y-fliij.,!  ai  -f-  *n  i  4-  y«{(f>7.i  -r  hi  )(*!  +  £)  -r  fta  +  &»  + 

+  4,S.i.iJ  —  «[2(Pr.i  +  P*.i)  +  +  +  +  01; 

c  ■»  2?3.i5ys  —  y*tf  -f  £».il  —  yfSa. i  —  2?6.i?)4'Ps.i! 

d  =■=  ?*&uy4  —  f  1 3n  +  fc,i )  ?y’  -r  i,;l r*6.i ?*  —  ^.iwnsa,  —  &ii?)  + 

~  >’!/»-[»,  a„52.i+  2 ^3.i(at  +  »«1  +  *i -  &•.>)]  +  ?  Pv)  + 

'  <*»sf  21  jjj.2  +  ijj)  —  4S;«.i  —  \ *th-  +  ®«Pm)  ~ 

*i«n?o.i  —  2laflir.i  —  ajp*.i)J; 

-  C2C,C4C,  -  C1CC,CV+  C|C?C,C»; 

V.  -C3C;CaC1,.~C4C\6;CI1 -r  C2C1C»Cl0+ClC,C,Cl#;  . 

32.5  —  -  C3CtC;C*  C,C3C4C»  —  C,CtC4C,; 

Pi. i  *■  C.C5CjCw;  “  — CtCjC(; 

Pi. i  ■“  CjCjCj)  -7-  C,CSC10  +  C,C,C,|  C|C}Cn; 

I'J  ~  ^iCjC,;  ^  C|C};  3?.j  =  CoCjCo CjCjCj -j- C|C*C»; 

3g.i  *=  CjC;C|  -7  CjCjC(i  -7-  C|CjC(;  j>s.j  =, — C|C|C|C^ 

?».?  *=  CtC3CQClt.  -f  C„(CjC4C,  —  C.CA  +  C.CA); 

P*./*=**lP*,y. 

• 

Characteristic  equation  (12)  has'  two  roots  a ^  and  0^. 
Accordingly,  the  general  solution  of  the  uniform  system  of  equations 
we  write  in  the  form: 

r-1  m.  n 

t-l  m.  n 

Tc  mm  Ji4—  —  —  (*»  I 

r  nz.-^+p; 

Pi  *“  ?-4  *  C|Cj,  j»j  =»  C|C]C|. 

P.i  P»  +  CjCjCjo-f*  CjCjCji  .  (13) 
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Expressions  of  analogous  type  are  written  for  w,.  The 

v 

particular  solution  of  heterogeneous  system  of  equations  (8)  we  look 
for  in  the  form: 


m,  « 


(/«  I.  n). 


Having  substituted  (14)  and  (7)  in  (8),  we  obtain  a  system  of 

heterogeneous  algebraic  equations  from  which  can  be  found  expressions 

for  coefficients  A*  and,  B*  through  coefficients  B  .  As  a  result 

c 

we  obtain  the  following  formulas : 

P?(  8), 

tf“R  eSto/"u,c')p:(»). 

zr  _ 

a--  (r,P, +2,Pl+zd«+2;?d,,+rl2,h+N)..v 

>,  -  CjCjCjC^,  P,  -  iC.C.C,. 

i1!  “  t|C|C|C(  -f  /CfQCsCn  -f  C|C,|).  (15) 

The  expression  for  E,  is  obtained  from  the  above  formula  for  F  with 

the  substitution  o  »  a  . 

c 

Thus,  the  general  solution  for  stream  function  is  written  in 
the  form: 

*,  -  R.  2J W*'-'  +  jf,  v>|ff.  ( 16 , 


Analogously  it  is  possible  to  give  the  solution  for  w.  Arbitrary 
constants  B1,  B 2  are  determined  according  to  initial  conditions. 
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The  process  of  their  determination  coincides  completely  with  that 
given  in  works  [2]  and  [3]  and  we  will  not  repeat  here  the  corre¬ 
sponding  conclusion. 

The  forecasting  scheme  described  above  was  computerized. 

During  calculations  values  of  coefficients  were  taken  as 
5  2 

fe2  *  5*10  m  /s,  fe  «  1^0  m,  and  the  remaining  parameters  were  the 
same  as  in  [2],  Forecasts  by  this  scheme  were  compared  with  fore¬ 
casts  according  to  the  scheme  given  in  [2],  As  a  quantitative 
estimate  of  the  quality  of  the  forecasts  we  used  the  relative 
forecast  error  (ratio  of  mean  absolute  error  of  forecast  to  mean 
absolute  variability).  Forecasts  of  pressure  at  sea  level  for  the 
North  Atlantic  were  analyzed  (zone  between  80°  west  longitude  and 
10°  east  longitude  and  between  20°  and  80°  north  latitude). 


Table  1. 


Original 

date 

i 

1 

1 

1  2 

E 

!  |  2 

3  III  1966 

1 

10.1 

1 

7.7 

I 

1.24 

096 

17  III 

9.4 

73 

138 

1.06 

•  22  III 

IO.fi 

«3 

132 

1,04 

27  III 

9.2 

73 

1,23 

0.92 

IS  IV 

l(U 

73 

1.44 

1.09 

18  IV 

113 

113 

1.13 

137 

221V 

133 

11.0 

1.61 

133 

27  1V  ! 

103 

Us 

1.40 

0.98 

291V 

73 

5.4 

.  0.96 

067 

«V 

i  11.1 

83 

1.23 

092 

13  V 

83 

.  7.2 

..  0.83  | 

!  072 

18  V 

33 

6,8 

132 

1.09 

aov 

93  ! 

1  -73 

1.12 

088 

23  V 

83 

7.7 

097 

091 

2S  V 

9.9 

7.7 

-  1.46 

1.16 

»v  ; 

M. 

73 

137  ; 

1  1.03 

1  VI 

83 

63 

132 

094 

Average ....  j 

\  '93 

73 

* 

136 

139 

Table  1  gives  estimates  of  the  quality  of  forecasts  for  four 
days  ahead:  6  -  mean  absolute  error  of  forecast  (decameter)  in  the 
considered  area,  E  -  relative  error.  The  magnitudes  of  the  estimate 
are  given  according  to  the  system  in  [2]  (column  1)  as  well  as 
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according  to  the  given  system  (column  2).  During  analysis  of  the 
data  of  Table  1  one  ought  to  consider  that  the  estimate  was  made 
for  regions  badly  illustrated  by  the  data  of  observations. 


The  data  of  Table  1,  and  also  immediate  comparison  of  forecast 
maps  composed  by  both  systems  testify  to  a  systematic  increase  in 
the  quality  of  forecast  using  the  new  system.  Prom  July  1966, 
ground  pressure  has  been  forecasted  using  the  tystem  described 
above. 
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FINITE-DIFFERENCE  ALGORITHM  OF  A  SOLUTION  OF  THE  VORTICITY 
EQUATION  FOR  THE  MIDDLE  TROPOSPHERE  OVER  THE 
NORTHERN  HEMISPHERE 


N.  V.  Isayev  and  M.  S.  Fuks-Rabinovich 


A  system  for  forecasting  the  geopotential  at 
the  middle  troposphere  using  a  space  (in  coordi¬ 
nates  xs  y)  finite-difference  approximation, 
which  provided  conservation  of  the  quadratic 
integral  features  is  considered  (i.e.,  vorticity, 
its  square  and  kinetic  energy);  integration  in 
time  follows  the  method  of  Adams.  Charts  of  the 
change  in  average  kinetic  energy  of  forecast 
fields  N^qq  in  time  are  constructed.  Results 

are  given  for  the  calculation  of  long-range 
forecasting  of  a  brief  length  of  time  before 
forecast  phenomenon  occurrence  according  to  the 
given  system. 


Introduction 


With  the  numerical  solution  of  forecast  equations  one  of  the 
greatest  difficulties  is  selecting  an  algorithm  which  would  provide 
stability  of  multistage  calculation.  In  view  of  the  complexity  of 
forecast  equations  and  the  large  volumes  of  original  and  inter¬ 
mediate  information,  it  was  impossible  to  show  earlier  a  method  of 
solving  the  forecast  problem  which  was  optimum  in  the  sense  of 
resistance  and  economy.  Itie  selection  of  such  a  solution  can  be 
made  only  on  the  basis  of  numerical  experiments  with  the  forecast 
model  being  investigated. 
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In  writing  a  finite-difference  model  of  linear  terms  one  must 
take  into  account  that  a  source  of  instability  can  be  inaccuracy 
in  the  approximations  of  both  the  space  and  time  derivatives.  In 
a  number  of  works  [8,  133  it  was  shown  that  a  certain  improvement 
of  the  quality  of  the  solution  can  be  obtained  by  using  a  higher 
approximation  (third  order)  of  the  first  and  second  space  derivatives. 
In  this  case,  however,  the  source  of  the  instability,  connected  with 
nonlinear  terms,  is  not  removed.  This  type  of  instability,  called 
nonlinear  instability,  is  shown  in  the  work  of  Phillips  [12].  The 
system,  having  linear  stability  (i.e.,  stability  investigated  on 
the  basis  of  linearized  equations),  was  far  from  always  suitable 
for  solving  nonlinear  tasks,  for  which  nonlinear  instability  can 
prevail. 


Approximation  of  Nonlinear  Terms  of  the 
Barotroplc  Vorticity  Equation 

An  attempt  to  remove  instability  connected  with  nonlinear 
terms  was  undertaken  in  [5].  It  was  proposed  to  use  a  finite- 
difference  system  of  calculating  the  Jacobian,  which  provided  the 
conservation  of  the  first  integrals  of  the  barotroplc  vorticity 
equation.  f 


Let  us  remember  that  for  the  barotroplc  atmosphere  the  absolute 
vorticity  of  velocity  fl  is  the  retained  quantity  and  the  vorticity 
equation  can  be  written  in  the  form 
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-0, 


(1) 


where  fl  •  +  1 1 

l  -  Coriolis  parameter 


u  and  v  -  horizontal  components  of  velocities 


> 


I 


at 


i+m£+9 


* 

*• 


On  the  other  hand,  the  general  kinetic  field  energy 


(2) 


does  not  change  with  time  under  the  condition  of  invariability  of 
values  of  the  stream  function  on  the  boundary  of  the  range  of 
integration  o.  In  integral  form  conditions  of  conservation  can  be 
written  in  the  form: 


£ —  Jj £(■*,  y,  /)do  —  const, 

• 

(3) 

5  — y.  t)d3  —  const, 

• 

(4) 

moreover,  the  constants  in  the  right  side  of  these  relationships 
are  a  function  of  time. 

It  is  obvious  that  during  the  transition  to  a  finite-difference 
type  of  forecast  equations  it  is  advisable  that  relationships  (3) 
and  (4)  be  executed  with  the  highest  possible  degree  of  accuracy. 
However,  the  majority  of  existing  difference  schemes  provides 
retention  of  only  part  of  the  integral  features  of  the  original 
equations  (for  example,  only  vorticity  velocity  ft  is  retained). 
Therefore  there  is  great  interest  in  the  work  of  Arakawa  [5],  which 
proposes  a  finite-difference  model  of  nonlinear  terms  of  the  vor¬ 
ticity  equation,  which  provided  an  integral  retention  of  the  vorticity 

velocity  fl,  square  of  the  vorticity  velocity  y.  /)  da  and 

» 

kinetic  energy  E,  i.e.,  all  integral  features  simultaneously. 

y 

In  the  work  of  Lilly  [8]  on  simplified  test  models  it  is  shown 
that  using  the  Arakawa  system  various  sources  of  nonlinear  instability 
are  mutually  compensated  and  nonlinear  instability  Is  exhausted. 

A  number  of  authors  [3,  9,  1*0  expresses  the  opinion  that  for 
purposes  of  short-range  forecasting  it  is  compulsory  not  to  use 
finite-difference  systems  in  which  the  Integral  features  of  the 
original  equations  are  accurately  retained,  but  It  is  possible  to 
use  a  system  where  retention  of  integral  features  Is  provided 
only  approximately.  But  for  long-range  forecasting  even  of  a  brief 
length  of  time  before  forecast  phenomenon  occurrence,  already  it  Is 
extremely  advisable  to  use  systems  which  provide  accurate  retention 
(at  least  within  the  framework  of  the  test  model)  of  integral 
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features  of  type  (3),  (4)  and  free  from  nonlinear  instability, 
because  during  multistage  calculation  in  long-range  forecasting 
nonlinear  instability  is  considerably  stronger  than  in  short-range 
forecasting. 

In  the  present  work  nonlinear  members  of  the  vorticity  equation 
were  represented  in  finite-difference  form  according  to  the  Arakawa 
system.  Following  [8],  let  us  describe  briefly  this  algorithm. 


Let  us  examine  the  vorticity  equation  in  the  solenoid  approach 
for  the  barotropic  atmosphere  (at  l  -  const) 


A,.  'A  -  —  —  4*- 

**  #»  #r  t* 


(5) 


where  c  *  Ms  ♦  -  stream  function,  J(c,  ♦)  -  Jacobian,  A  -  two- 
dimensional  Laplacean  operator. 


We  introduce  the  designations: 


«-  ~ (/■*! x  -r  4p  |  -  Mjr--y-)|. 

Jf » «  t\  v  |>U  -r  i  jt)  -  fix  -  a.rjJ , 


where  P(x)  is  any  of  the  considered  functions. 

In  (6)  let  us  write  down  three  such  difference  expressions 
for  Jacobian  J: 


A  -  -yu.** )*  +  M- (7) 


According  to  [8],  each  of  these  relationships  provides  vorticit 
retention,  the  vorticity  square  and  kinetic  energy  respectively. 
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If,  however,  we  put  together 


Ja**  Jf¥  Jt)-  ( 8 ) 

then  this  finite-difference  model  will  guarantee  retention  of  all 
three  shown  features  simultaneously. 

Being  limited  in  the  approximation  of  nonlinear  terms  to  the 
second  order  of  accuracy,  we  have  the  following  working  formulas 
for  determination  of  J,  (see  Pig.  1): 

A 

J  'I*'*  —  *iK —  — 

-<v.  -  '..Xy,  - 

~  J’*-  V:)  —  •.•4*'.  \)  4"  *<(v»  —  + 

(9) 

where  h  -  grid  interval. 


Pig.  1.  Grid- 
template  for 
calculating  non¬ 
linear  terms  (h 
grid  Interval). 


Approximation  of  Derivative  in  Time 


In  forecasting  systems  and  in  models  of  the  overall  circulation 
are  applied  various  methods  of  approximating  the  derivative  in  time. 
Along  with  the  method  of  central  differences,  various  authors  use 
the  systems  of  Lax-Wendroff  [7],  Mlyakoda  [11],  Matsuno  [10],  Hoyne 
[6],  and  others.  According  to  test  calculations  In  a  work  of  Lilly 
[8]  a  good  index  of  stability  goes  with  the  method  of  Adams  [1]. 

Per  integration  In  time  this  work  tested  the  system  of  Adams, 
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Miyakoda,  Hoyne,  Matsuno  and  the  method  of  central  differences. 

The  results  which  were  best  and  very  close  In  quality  were  obtained 
with  the  system  of  Adams  and  Miyakoda. 

In  Table  1  are  estimates  of  various  methods  of  integration  in 
time.  Along  with  the  change  in  kinetic  energy  for  the  various 
systems  the  mean  relative  forecast  errors  by  these  systems  are  given. 
Averaging  took  in  75  points.  Furthermore,  the  amounts  of  divergence 
between  the  different  systems  are  given.  The  best  results  were  with 
the  systems  of  Adams  and  Miyakoda  (however,  the  divergences  between 
them  are  small).  Both  these  systems  are  more  accurate  than  the 
method  of  central  differences  (and  the  methods  of  Matsuno  and  Hoyne). 


Table  1.  Comparison  of  certain  methods  of 
integration  in  time. 


Method 

-  -  V“~ 

Formula 

7-1 

7—2 

7—3 

ir»|  ^ 

1 

\r\ 

'l-lr,i|ri 

■c. 

! 

j 

Miyakoda 

r  i. t  o  . 

!.l  |  A7I 

f 

1 

7.3 

U 

xs 

(Mil 

M 

3.7 

0.90 

12.9 

Adams 

1  \ 

:  *-i 

;  ^  [ 

1  •  *l 

k 

i 

iron 

7.2 

*  • 

2.*  O.W 

10.1 

W 

0L« 

:x* 

■j 

Jr'**"*,-**  -  /\V 

1 

Hoyne 

i.v  fi.aa 

I 

K.2 

5.1 

a96 

\12 

to 

1.06 

14  9 

Central 

differ¬ 

ence 

x*-'  -ifM 

1 

■ 

2.1  VS 

V.4 

I.W 

14.2 

?t> 

I.2u 

■M 

Note.  xn  -  value  of  geopotential  on  the  n-th 
time  interval  /*-  |^j  -  value  of  deriva¬ 
tive  in  time  from  geopotential  on  the 
n-th  interval;  |f/j  -  absolute  deviation 
in  the  kinetic  energy  in  percent,  e  - 

mean  relative  forecast  error;  !a| 

’  cp 

mean  absolute  error  of  forecast  in  deca¬ 
meters;  At  -  time  interval;  T  -  number 
of  days;  8  -  parameter  in  system  of 
Miyakoda, 
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In  using  the  -methods  of  Adams  and  Miyakoda  kinetic  energy  is 
retained  to  a  greater  extent  in  the  course  of  calculation  than 
otherwise.  The  method  of  Adams  is  distinguished  also  by  its  sim¬ 
plicity,  it  is  handy  for  using  on  a  computer  and  requires  keeping 
only  one  additional  field  (as  the  method  of  central  differences).  In 
view  of  the  greatest  effectiveness  and  economy  we  stayed  with  the 
method  o-f  Adams  and  used  it  in  this  work  for  approximation  of  the 
derivative  in  time. 

We  will  give  the  derivation  of  the  Adams  extrapolation  formula, 
limiting  ourselves  to  an  approximation  of  the  second  order  of  accuracy. 

Decompose  functions  and  ipn  (values  of  stream  function  on 

the  n  +  1  and  n  time  interval  respectively)  in  a  Taylor  series 
in  At  and  2At 


•V'  -  -t  -  V  +  o(A 
V»- 1  “  V  I  +  +  0(A 


(10) 


We  multiply  the  first  and  second  relationships  from  (10)  by 
certain  real  numbers  a  and  0  respectively  and  combine  the  results 


(*'t  SOi'/i-m  -  i  -f  i**  +  2? )Af  + 


(11) 


Considering  the  third  relationship  from  (10),  let  us  rewrite 
formula  (11)  in  the  form: 


(*  +  1  =(Jr  +  (a  f,  +  ?  + 


(11 '  ) 
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In  order  to  obtain  an  approximation  of  the  second  order  of  accuracy, 

p 

the  term  of  relationship  (11*)  containing  A t  must  be  equivalent  to 
zero.  For  this  condition  it  is  necessary  that 


whence 


Substituting  these  values  into  expression  (11'),  we  obtain  the 
Adams  formula 

V '  -  *.  ■ +  (t  V.  -  T  Cr) * t + (12) 


Parameters  of  the  Forecast  System 

On  the  basis  of  the  algorithm  given  in  previous  sections  a 
system  was  constructed  to  forecast  the  geopotential  for  00  over 
almost  all  the  northern  hemisphere  (Fig.  2).  As  in  [4]  the 
barotropic  vorticity  equation  in  the  quasi-geostrophic  approach  is 
written  in  the  following  form: 


(— j-  A//+  /,  H\ 


(13) 


where  H  -  geopotential,  l  =  l($)  -  Coriolis  parameter;  m  =  m(<t>)  - 
parameter  of  the  increase  in  stereographic  projection;  $  -  geographi¬ 
cal  latitude, 


(a.  *)« 


•  hi  dt  fed* 

W  Ty  15j~~5x 


(14) 


The  original  geopotential  field  H  is  given  in  2181  mesh  points 
(Fig.  2),  having  the  form  of  an  octagon  (51  points  in  the  longest 
columns  and  lines  of  the  octagon  and  23  points  in  the  shortest) 
with  the  center  of  symmetry  at  the  north  pole.  Nonlinear  terms 
were  approximated  according  to  the  Arakawa  system,  and  the  derivative 
in  time  by  the  method  of  Adams.  Equation  (13)  relative  to  3 H/dt 
was  solved  by  the  extrapolation  method  of  Liebmann.  In  this  case 
at  the  edge  of  the  octagonal  area  a  condition  of  invariability  of 


17 
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Fig.  2.  AT^00  forecast  map  for  0300, 

24  December  1963  according  to  data  for  0300, 

22  December  1963.  The  solid  line  outlines  the 
octagon  along  which  the  forecast  was  resolved. 


geopotential  with  time  was  assumed.  During  the  calculation  no 
smoothing  operators  were  applied.  The  interval  in  coordinates 
was  390  km  at  a  latitude  4>  ■  60°,  the  time  interval  At  =  3  hours. 
This  system  resolved  forecasts  of  the  H^qq  fields  for  periods  of 
24,  48  and  72  hours.  Analysis  of  the  forecast  fields  will  be 
conducted  below. 

Calculation  and  Adjustment  of  the  Energy 
of  Forecast  Fields 


Retention  of  the  integral  features  provided  by  the  Arakawa 
system  nevertheless  is  distorted  during  the  calculation  of  a  fore¬ 
cast  for  a  long  period.  The  reason  for  this  is  not  only  the  errors 
of  extrapolation  in  time,  but  a  whole  series  of  circumstances. 

It  is  widely  known  that  basic  factors  limiting  the  earliness  of 
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the  forecast  model  are  first  of  all  the  simplifications  accepted 
in  the  derivation  of  forecast  equations  (such  as,  for  example, 
adiabaticity ,  quasi-geostrcphicity ,  and  others),  and  also  inaccuracies 
in  the  initial  information  and  boundary  conditions.  (In  forecasting 
for  the  hemisphere  or  sphere  specific  difficulties  appear  during  the 
examination  of  phenomena  in  the  equatorial  range.  Furthermore, 
limited  computer  capabilities  force  us  to  resort  to  economical 
storage  of  intermediate  information  during  the  calculation,  which 
leads  to  a  rounding  of  the  intermediate  values  and  limits  the 
accuracy  of  calculations  because  of  the  accretion  of  rounding 
errors . ) 

Here  we  will  describe  a  numerical  experiment  conducted  to 
determine  the  interval  of  earliness  of  forecast  of  the  given  system. 
Basic  criteria  for  determination  of  forecast  successes  were  of  course 
synoptic  analysis  and  quantitative  estimation;  however,  in  addition 
the  following  was  done.  . 

On  every  interval  in  time  during  the  calculation  of  the  fore¬ 
cast  the  quantity 


was  figured,  where  En,  EQ  are  values  of  the  kinetic  energy  fields 
calculated  by  formula  (*0  on  the  n-th  and  zero  time  intervals.  The 
quantity  W  shows  how  much  (In  percent)  the  kinetic  energy  of  the 
field  on  the  n-th  time  interval  changed  in  comparison  with  the 
initial. 

In  vorticity  equation  (13)  a  term  considering  horizontal 
turbulent  exchange  was  introduced,  whereupon  this  equation  took 
the  form 


ATF-(-TrA//+/* 


(16) 


where  v«=^+7-d~  is  a  plane  of  Hamiltonian  operator;  i,  J  -  unique 
vectors . 
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With  the  aid  of  parameter  v  the  kinetic  energy  of  the.  forecast 
fields  was  adjusted,  namely,  in  the  process  of  calculating  the 
forecast  the  quantity  v  was  selected  such  that  the  deviations  in 
energy  W  would  be  located  within  certain  limits 


(17) 


where  6  was  assigned  as  0.5J. 

Until  inequality  (17)  was  disturbed,  during  the  calculation  of 
the  forecast  the  value  v  ■  0  was  used.  When  inequality  (17)  was 
disturbed,  v  was  determined  from  the  relationship 

(18) 

where  k  -  an  empirically  selected  coefficient,  fe  -  number  of  interval 
on  which  occurred  the  disturbance  of  (17).  If  on  subsequent  intervals 

(17)  held,  then  v  did  not  change.  Such  intermittent  inclusion  into 
calculation  of  control  parameter  v  allows  regulating  the  kinetic 
energy  of  the  obtained  fields  and  achieve  its  change  within  the 
limits  of  2%  (for  a  lower  value  of  6  in  (18)  it  is  possible  to  lower 
this  limit  somewhat)  with  a  practically  unlimited  number  of  time 
intervals.  The  change  calculated  by  actual  data  for  three  days 
likewise  is  located  within  the  limits  of  several  percent. 

The  authors  calculated  the  geopotential  field  for  20  days; 

the  chart  of  the  change  W(t)  in  time,  obtained  as  a  result  of  numerical 
experiments  using  control  parameter  v,  which  was  determined  by  formula 

(18)  upon  disturbance  of  condition  (17),  is  given  In  Pig.  3b. 

Changes  In  values  of  v  during  the  numerical  experiment  are  shown 
in  Pig.  3a.  We  see  that  up  to  five  days  the  forecast  Is  resolved 

2 

at  very  small  values  0  <  W  £  2 ’10  m  /s ,  and  in  the  subsequent 

period  the  utilized  values  of  v  fluctuate  from  zero  ano  approximately 
5  2 

to  5*10  m  /s,  which  agrees  well  with  the  values  of  v  used  usually 
in  calculations  (see  for  example,  [2]). 

The  most  interesting  feature  of  the  chart  of  V(t)  (Pig.  3b) 
is  the  fact  that  after  five  days  the  quantity  |l/|  somewhat  increases 
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Pig.  3.  Charts  of  the  change  in  parameter 
lg  (y  +  10)  in  time  (a)  and  the  quantity  En/E Q 

in  time  (b).  E  r-  field  energy  on  the  n-th  time 
n 

interval.  -  energy  of  initial  field. 

/ 

and  the  change  of  V(t)  becomes  periodic.  This  periodicity  Indicates 
that  at  approximately  the  fifth  day  a  large  role  begins  to  be  played 
by  the  balance  of  purely  stray  perturbations  and  damping  factors, 
realizable  with  the  help  of  control  parameter  v.  Visual  and  quanti¬ 
tative  estimation  of  forecast  fields  likewise  shows  that  four  to 
five  days  is  the  limit  period  for  good  Justification  of  a  forecast 
based  on  this  system. 

Such  a  mechanism  of  estimating  the  earliness  of  forecast 
apparently  may  be  applied  to  other  forecasting  systems. 

Analysis  and  Estimation  of  Prediction  Fields 

As  already  was  indicated  above,  the  calculation  of  forecasts 
was  conducted  for  v  ■  0  (in  contrast  with  numerical  experiments  on 
determining  the  interval  of  earliness  of  forecast,  examined  in  th_ 
previous  section). 


During  the  calculation  of  forecasts  simultaneously  by  formula 
(1*0  the  quantity  W(t)  was  calculated  in  order  to  be  able  to  observe 
the  degree  of  retention  of  kinetic  energy  in  the  fields  being 
forecasted.  For  a  period  of  two  to  three  days  the  quantity  W(t) 
was  within  limits  of  2-3. 5%,  which  indicates  good  retention  of  the 
integral  features  for  such  a  period  of  calculation.  As  a  comparison 
let  us  say  that  when  using  the  algorithm  of  forecast  calculation 
given  in  [4],  the  quantity  V(t)  changed  after  three  days  by  15-20J, 
which  indicates  nonretention  of  the  integral  features  as  well  as 
a  large  degree  of  smoothing  of  the  forecast  fields. 

As  an  example  let  us  give  the  results  of  forecasting  for 

48  hours,  calculated  according  to  original  data  for  0300,  22  December 
1963  (Fig.  4).  The  actual  and  forecast  maps  for  0300,  24  December  1963 
are  given  in  Figs.  2  and  5.  (Let  us  note  that  we  selected  the  most 
difficult  example  during  the  period  of  the  two-year  test  in  a  system 
based  on  limited  territory.) 


Pig.  5.  AT500  actual  "^P  for  0300,  24  December 

1963. 


We  see  that  values  of  the  geopotentlai  in  centers  were  predicted 
rather  well,  although  a  certain  lag  from  their  actual  movement  is 
observed.  Cyclogenesis  over  the  Pyrenees  and  the  formation  of  a 
ridge  over  western  Germany  (where  there  is  a  trough  on  the  original 
map)  were  reflected  in  the  forecast.  The  absence  of  smoothing  and 
retention  of  integral  features  guarantees  that  the  forecast  maps 
will  have  pressure  gradients  close  to  the  actual  ones. 

Quantitative  estimates  of  the  example  are  given  in  Table  2. 

Over  the  illustrated  regions  relative  error  e  »  O.83  (over  Europe 
c  *  0.72),  whereas  forecast  estimates  calculated  by  the  simpler 
system  described  in  [4]  are  noticeably  worse:  c  «  1.2-1. 4.  (In 
Table  2  estimates  were  resolved  over  all  points  of  the  field.) 

An  estimate  of  calculation  results  for  a  series  of  forecasting 
examples  from  this  system  are  given  in  Table  3.  The  estimate  was 
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Table  2.  Estimate  of  forecast  success  for  48  hours 
according  to  original  data  for  0300,  22  December  1963 
(I  -  finite-difference  model  in  work  [4];  II  -  in  this 
work ) . 
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Table  3*  Estimates  of  forecasts  for  24,  48  and  72  hours 
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over  75  points  (£Cp)>  as  far  as  possible  evenly  distributed  through 
the  forecast  range,  and  also  through  12  points  (c^),  located  in 
the  North  Atlantic  and  on  the  European  territory  of  the  USSR. 

One  ought  to  note  that  the  examples  of  forecasts  using  raw 
data  for  0300  15-19  December  1965  were  specially  selected  for  the 
purpose  of  testing  the  system  on  complex  synoptic  situations. 

Analysis  of  a  series  of  examples  of  forecasting  (12  cases), 
calculated  by  this  system,  allows  us  to  make  the  following  quali¬ 
tative  conclusions.  First  of  all  one  ought  to  note  that  such 
Important  phenomena  in  reconstruction  of  fields  as  cyclo-  and  anti¬ 
cyclogenesis  frequently  are  provided  in  forecasts  (cases  of  cyclo- 
and  anticyclogenesis,  conditioned  basically  by  dynamic  factors). 
Furthermore,  the  change  in  pressure  in  the  centers  of  the  baric 
formations,  and  also  values  of  gradients  in  the  pressure  field  are 
correctly  predicted.  At  the  same  time,  a  certain  (of  the  order  of 
20t)  error  in  predicting  movements  of  the  centers  of  the  baric 
formations  is  observed.  (The  last  drawback  can  be  in  some  measure 
removed  by  increasing  the  order  of  approximation  of  the  nonlinear 
terms. ) 


In  conclusion  let  us  note  that  this  finite-difference  algorithm 
can  be  applied  even  to  baroclinic  forecast  systems. 
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NUMERICAL  SOLUTION  OP  THE  BALANCE  EQUATION 
IN  THE  FRAMEWORK  OP  A  QUASI-SOLENOIDAL 
FORECAST  DIAGRAM  OP  THE  OEOPOTENTIAL 
FOR  THE  NORTHERN  HEMISPHERE 


I.  G.  Sitnlkov  and  S.  0.  Krichak 


A  finite-difference  system  of  solving  the 
balance  equation  for  the  northern  hemisphere  is 
given.  A  series  of  characteristic  features 
found  with  the  solution  of  this  equation  for 
a  belt  of  low  latitudes  is  shown.  The  obtained 
values  of  the  stream  function  are  used  to  fore¬ 
cast  the  geopotential  field  at  the  middle  level 
of  the  troposphere  for  72  hours.  An  example 
of  the  forecast  is  presented. 


Various  questions  connected  with  the  technology  of  numerical 
solution  of  the  balance  equations  have  been  illustrated  in  a  number 
of  works  (for  example,  [1,  3,  5,  6,  7,  8,  9]).  Most  frequently  in 
this  case  the  problem  is  solved  for  a  comparatively  small  area.  Of 
the  mentioned  work.’,  only  In  [8]  is  a  method  stated  for  solving  this 
equation  relative  to  the  northern  hemisphere.  'this  work  establishes 
a  number  of  characteristic  features  found  with  the  solution  of  the 
balance  equation  within  the  framework  of  a  quasi-solenoidal  system 
of  forecasting  the  geopotential  on  the  northern  hemisphere. 


1.  The  balance  equation  recorded  in  a  local  isobaric  system 
of  coordinates  fx,  y,  p)  for  a  polar  stereographic  projection  has 
the  form  [9] 
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2m»  U  *4,  ,»  **  #41. 

m  *  ■  T[( dxdy )  a.v-  j +  /  i  a* 


a/  .  a^-  a/\ 

?Fr 


^a//. 


Here  tJ»  -  stream  function,  -  height  of  isobarie  surface,  m  - 
"parameter  of  distortion"  of  stereographic  projection,  Z  -  Coriolis 
parameter,  g  -  acceleration  of  gravity. 


It  is  known  that  equation  (1),  considered  as  the  equation  for 
determination  of  values  of  4s  will  belong  to  the  family  of  Monge- 
Ampere  equations  referring  either  to  elliptical  or  hyperbolic  type 
depending  on  the  relationship  of  coefficients  [*)].  A  number  of 
authors  [8,  9],  making  use  of  the  so-called  Petersen  conversion, 
reduce  equation  (1)  to  the  form 


/Jl3A  i  e=  /  l-  2 g  A  H  -f-  /H*(  - 


djfi 


4? 


+ 


(2) 


To  solve  equation  (2)  usually  the  method  of  series  approxima¬ 
tions  is  used.  In  this  case  equation  (2)  is  considered  as  a  Poisson 
equation  relative  to  ip,  considering  that  values  of  ip  in  the  right 
side  are  known  and  are  determined,  for 'example,  from  the  previous 
approximation. 


We  will  use  as  the  initial  approximation  for  the  solution  of 
equation  (2)  the  relationship 

po-f//,  (3) 

where  l  -  the  value  of  l  on  any  fixed  latitude  <pc  (we  took  =  45°). 
Then,  substituting  (3)  in  equation  (2),  we  can  after  solving  (2) 
determine  according  to  4>q  a  new  approximation  ipj ;  substituting  it 
In  the  right  side  of  (2),  we  obtain  ipp ,  etc. 

As  a  boundary  condition,  following  Miyakoda  [8],  we  take 
relationship  irp  =  '—tflrp.  where  again  l-/(tpo),  ?>o*»450. 
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The  obtained  field  ip  was  used  as  the  initial  field  for  the 
quasi-solenoidal  forecast  in  terms  of  a  barotropic  model  on  the 
northern  hemisphere  (system  proposed  in  [2]). 

The  forecasts  of  were  for  24,  48  and  72  hours.  At  the  end 
of  each  period  the  field  was  located  by  means  of  so-called 

"inversion"  of  the  balance  equation,  i.e.,  of  the  solution  relative 
to  3  for  known  distribution  of  ip. 

2.  For  a  solution  of  equation  (1)  we  must  require  that  this 
equation  be  elliptical  in  all  points  of  the  grid.  As  experiment 
shows,  the  criterion  of  ellipticity  of  equation  (1)  having  the  form 

r  -  h+  n  -  2(4t  £ + %  -£•)  >o,  (,) 

is  observed  for  synoptic  objects,  at  least  for  moderate  latitude, 
and  with  a  grid  interval  of  no  less  than  250  km  [4], 

Usually  the  last  term  of  (4)  in  absolute  value  is  considerably 
less  than  the  first  two  terms,  therefore  frequently  (see,  for 
example,  [8,  9])  a  simpler  criterion  of  ellipticity  is  used. 

r  =  2zm!LH4-  /-\>0.  (5) 

This  work  assumed  criterion  (5).  On  Fig.  1  is  shown  the 
distribution  over  the  northern  hemisphere  of  points  with  r  <  0  for 
one  of  the  analyzed  examples.  The  same  figure  represents  the  E 
field.  During  the  calculation  of  T  values  of  Laplacian  A H  were  taken 
with  fie  =  390  km  at  a  latitude  of  60°  (at  a  latitude  of  10°,  to  which 
the  extreme  grid  points  reached,  the  grid  interval  was  about  200  km). 
The  number  of  grid  points  in  which  equation  (1)  is  hyperbolic  attains 
here  400,  i.e.,  18J?  of  all  grid  points.  In  this  case  the  over¬ 
whelming  majority  of  these  points  falls  in  low  latitudes,  <(>  <  30°. 

For  the  belt  of  latitudes  <f>  >  45°  the  number  of  points  of  the 
hyperbolicity  of  equation  (1)  is  25. 
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Fig.  1.  Distribution  of  points  of  hyperbolicity 
of  equation  (1)  for  the  northern  hemisphere  for 
a  square  grid  with  an  interval  6s  «  390  km  at 
latitude  <j>  *  60°  and  j?^00  field  at  0300  17  Septem¬ 
ber  1965.  j 

As  one  would  expect,  a  considerable  number  of  points  for  which 
T  <  0  is  connected  with  the  central  ranges  of  deep  anticyclones 
(mainly,  subtropical),  where  AJ7  Laplacians  take  large  negative 
values . 

It  is  possible  to  fulfill  the  criterion  of  ellipticity  of  (5) 
for  all  grid  points  by  correcting  the  values  of  bH  in  those  areas 
where  r  <  0.  This  was  done  in  the  following  manner.  Assume  that 
at  point  0  (Fig.  2)  r  ■  Tg  <0.  Write  criterion  of  ellipticity  (5) 
in  finite-difference  form 

(6) 

(here  SH ,  for  example,  at  point  0  is  expressed  as  (6H) Q  »  H 1  +  H 2  + 

+  H  -  +  H A  —  4H n)  . 

3  4  0 
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Pig.  2.  Finite- 
difference  approxima¬ 
tion  of  first  and  second 
derivatives  in  neighbor¬ 
hood  of  point  0. 


Divide  both  sides  of  (6)  by  m2/(6a)2.  We  obtain 

^*5?-  Wf  f  >  0:  (6a) 

Let  now  yq  «  0  be  the  new  value  of  y  at  point  0.  Add  to  each 
of  the  values  of  y .  (i  *  1;  2;  3;  *0  the  quantity  yn/4  and  determine 
new  values  of  y.  at  points  1-4  as  y.  *  y.  +  (y,,/4).  This  means 

t  Z  Z  U 

that  the  mean  value  of  y  in  nearest  neighborhood  of  point  0  becomes 
constant. 

We  consider  now  that  the  values  of  field  B  at  points  0,  1,  2, 

3,  4  are  unknown  (designate  them  etc.),  and  determine  them  from 

the  following  system  (see  Pig.  2): 

H,  +  H,  +  Ht+Ht-  AH.  -  -^(t, 

H„  +  H%  -f  H,  +  H'  -  AH,  - -£-(7,  -  L, ); 

«o  +  //» •+  //,«  +  «.  -  4//,  -  -  Lt); 

/>,  +  //,  +  H„  +  //t  -  4//,  -  -jL-ta  _  i  j. 

Mi  "f  W,  -f  H,  -f  H'  —  AH, »  ^«)» 

Solving  system  (7),  we  find  new  values  of  HQt 
such  that  the  value  of  y  at  the  investigated  point  (point  0)  is 
equal  to  zero.  The  described  method  is  close  to  the  method  used 
by  Shuman  [9], 
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We  agree  to  call  numerical  values  of  y^  in  (6a)  the 
"hyperbolicity "  at  an  arbitrary  field  point  i.  Experiment  showed 
that  by  the  examined  means  it  is  possible  to  blur  the  hyperbolicity 
at  the  given  point  over  the  nearest  neighborhoods  of  this  point. 
Having  applied  the  described  procedure  to  all  points  of  the  grid 
where  y.  <  0  (except  the  boundary),  we  decrease  the  absolute  value 
of  both  maximum  hyperbolicity  ym_u„  and  average  hyperbolicity  y 

naKC 

over  the  whole  field  (  where  the  summing  is  over  all  B 

Inside  field  points  in  which  y  <  0)  because  of  the  more  precise 
determination  of  values  of  H  at  points  where  earlier  it  was  y  2:  0. 

For  a  further  decrease  in  the  values  of  y  and  y  it  is  necessary 

nOnC  Cp 

to  make  several  circuits  of  the  entire  map. 

We  considered  the  task  completed  if  | YMaKC I  <  0*25*  Note  that 

after  the  first  circuit  of  the  map  tne  quantity  I >'MaK c  I  was  usually 

10.0-30.0.  For  field  B  shown  in  Fig.  1  after  the  first  circuit 

YMaHC  =  -20-3»  Ycp  *  -2.75;  after  the  tenth  circuit  yMaKC  =  -0.23; 

y  *  -0.02.  In  this  case  the  maximum  difference  between  the  values 
cp 

of  original  and  corrected  field  H  did  not  exceed  3  decameters  in 
even  one  point .  The  greatest  corrections  were  necessary  on  low 
latitudes,  and  regions  with  weak  meteorological  treatment. 

The  modified  field  H  was  used  as  the  initial  field  for  solving 
the  balance  equation  relative  to  ^ 

3.  Let  us  write  equation  (2)  in  finite  differences  for  point  0 

2  2 

(Fig.  2),  having  multiplied  the  left  and  right  sides  by  (fie)  / mQ 
(m0  -  value  of  "distortion  parameter"  m  at  point  0).  We  have 

V +  t*  +  +  ^ 

*  -f  +  +  4//,)+  -!(*»+  fc- 

*4)*-  ^-Kti -*.><*»  -  /,)  +  (♦<  -  *)<*,-/,))  (8) 
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2 

(here  As  -  value  of  step  of  grid  in  decameters:  As  =  (£a).‘10  ,  where 
6e  *  390  km;  ZQ  -  value  of  Coriolis  parameter  at  point  0).  During 
the  calculation  of  finite-difference  models  of  the  second  derivatives 

-5^-,  first  differences  were  used  with  an  interval  /26a 

over  the  grid  (x^,  y ^),  rotated  45°  in  comparison  with  the  original. 
Thus,  at  point  0  we  have  (Pig.  2): 


l*i  *yi 


.JLAeA 

ly,  »  ft* 

ijr,  TFiT 


:±±±~*b  . 

»»)• 

,k±iiz2fe.. 

*»*>» 


(9) 


Such  a  recording  of  the  second  derivatives  is  close  to  "method  C" 
analyzed  in  [8].  It  is  the  finite-difference  expressions  of  first 

derivatives  of  type  *j£  (q  •  x,  y)  taken  with  an  interval  of  26a, 

Laplacian  A E  with  an  interval  of  /26a,  Laplacian  Ai|>  in  the  left 
side  of  (8)  with  an  interval  of  6s. 


For  calculations  from  formula  (8),  furthermore.  It  was  con- 

_ji 

venient  to  multiply  both  sides  of  the  equation  by  10  ,  considering 

-4  2 

new  function  if»*  *  ^*10  .  The  value  of  function  In  dam/s  ,  as 

#•10““*  - 

is  simple  to  see,  for  example,  from  expression  — j — H  (g  *  0.98 

p 

dam/s  ),  has  the  same  order  of  magnitude  as  the  values  of  H  in 
decameters . 


To  solve  system  of  equations  (8)  in  all  inside  points  of  the 

grid  the  extrapolated  Liebmann  method  was  applied,  modified  in 

1  2 

the  following  manner.  Namely,  the  series  approximations  , 

•j 

ip  ,  etc.,  developing  in  the  solution  process,  were  substituted  in 
the  right  sides  of  equations  (8),  not  every  time,  but  only  after 
the  fulfillment  of  some  fixed  number  of  iterations  N.  We  took 
N  *  10.  Thus  new  fields  F  were  formed  after  the  10,  20,  30th,  etc., 
iteration,  and  the  entire  process  of  calculation  consisted  of  two 
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cycles  "imbedded'’  in  one  another.  Calculations  were  conducted 
up  to  the  achievement  of  the  final  accuracy  of  6,  in  other  words, 

it  was  necessary  that  | -ts- 

This  technology  of  solving  system  of  equations  (8)  is  close  to 
the  "accelerated"  method  (fast  method),  given  in  work  [8,  9].  For 
an  accuracy  of  i  ■  t  dam  it  was  necessary,  as  a  rule,  to  have  80-100 
iterations  (the  field  of  right  sides  was  calculated  only  8-10  times). 
This  calculation  occupied  about  30  minutes  of  machine  time. 

In  the  course  of  a  solution  in  all  grid  points  where  the  sub¬ 
radical  expression  in  formulas  (8)  is  less  than  zero  (observance 
of  criterion  (5)  leaves  this  possibility),  it  is  equated  to  zero. 

It  was  clarified  that  the  number  of  points  in  which  it  is  necessary 
to  produce  this  correction  is,  as  a  rule,  150-200  (7-9%  of  all  grid 
points),  moreover  these  points  are  distributed  mainly  on  belt  of 
low  latitudes.  This  effect  can  be  connected  with  the  fact  that  in 

low  latitudes  the  term  in  the  subradical  expression 

In  the  right  side  of  equation  (2)  sometimes  acquires  large  positive 
values.  It  is  possible  that  a  correction  of  field  B  in  accordance 
with  criterion  of  ellipticity  (4),  instead  of  (5),  would  bring  down 
the  number  of  points  in  which  an  artificial  change  In  the  size  of 
the  subradical  expression  is  required. 

We  must  note  one  circumstance  leading  to  a  reduction  of  the 
rate  of  convergence  during  calculation  for  low  latitudes.  As 
experiment  showed,  the  calculated  values  of  <f/*,  numerically  equal 
in  the  geostrophlc  approach  to  values  B  dam,  in  the  quasi-solenoidal 
approach  seem,  as  a  rule,  more  than  B.  This  one  can  be  seen  well 
in  Fig.  3.  In  moderate  and  high  latitudes  */•  gradients  are  less 
than  B  gradients,  but  in  low  latitudes,  under  condition  (3),  fixing 
values  cf  \ p*  on  the  boundary,  gradients  grow  in  comparison  with 
H  gradients.  Such  a  configuration  of  field  in  low  latitudes 
leads  to  a  drop  in  the  rates  of  convergence  of  the  iteration  process 
(in  comparison  with  the  initial,  geostrophlc,  approach). 
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Fig.  3.  Cut  along  45°  west  longitude  and  135° 
east  longitude^for  the  ^ield  at  °300 

17  September  1965  and  the  4>§  field,  obtained 
by  solving  the  balance  equation. 


4.  Let  us  stop  at  the  question  about  the  inversion  of  the 
balance  equation.  Inasmuch  as  equation  (1)  relative  to  function 
H  is  Poisson  equation.  Its  numerical  solution  does  not  represent 

p 

difficulty.  Multiply  equation  (1)  by  l/m  g,  exchange  the  left  and 
right  sides  and  write  It  in  finite-difference  form  (see  Pig.  2) 

+  +  +  Ht  -  4H,  — (r-(«  +  f.  +  *  +  * - - 

U 

-  ♦*, -  w- -  ( f. +  *?■ - 20(15* +  « - *01  +■• 

+  *rr  U  -  b)(l  -  l)  +  (tf  -  ♦$/.  -  OL  ( 1 0 } 

Ug 

■'*'  "•  p  «•  ij  ^ 

Here  l  »  Z'lO  ;  g  *  0.98  dam/s  ;  ij/*  -  ^*10  ;  expressions  .  ggg^ 

are  solved  according  to  formulas  (9).  A  system  of  equations  of  form 
(10)  we  solve  by  the  extrapolated  Liebmann  method.  To  achieve 
an  accuracy  of  6Q  s  0.5  dam,  as  a  rule  about  10  iterations  are 
necessary.  Inversion  of  the  balance  equation  after  a  forecast  for 
2*1,  48  and  12  hours  occupies  in  overall  complexity  about  10  minutes 
of  machine  time. 
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Inversion  of  the  balance  equation  was  used  also  in  the  next 
experiment.  A  system  of  equations  of  form  (8)  was  solved  relative 
to  ij>#,  and  in  terms  of  the  found  values  the  inversion  immediately 
was  executed  according  to  formulas  (10).  The  inverted  field  must, 
generally  speaking,  completely  repeat  the  original  field  B  (except 
those  points  in  which  either  the  criterion  of  ellipticity  (5)  was 
not  executed,  or  subradical  expressions  in  (8)  were  equated  to  zero). 
The  absence  of  iteration  indicates  either  deficiencies  of  the 
finite-difference  recording  of  (8)  or  (10),  or  inaccuracy  in  the 
direct  or  inverse  solution  of  the  balance  equation. 

Comparison  of  the  profiles  of  original  and  inverted  B  fields 
shows  that  the  agreement  is  entirely  satisfactory. 

5.  A  finite-difference  system  of  barotropic  prediction  of 
the  geopotential  on  the  northern  hemisphere  in  the  quasi-geostrophic 
approach  is  described  in  [2].  This  system  was  used  by  the  authors 
to  forecast  the  values  of  stream  function  4/*,  obtained  as  a  result 
of  solving  the  balance  equation.  The  balance  equation  was  inverted 
following  formula  (10).  In  Pig.  4  are  given  the  contours  of  the 
area  for  which  the  forecasts  were  ir.ad4  ^y  the  system  (grid  Interval 
6b  *  390  km  at  latitude  60°),  and  the  forecast  field  for  1500 

hours,  10  October  1961,  calculated  48  hours  in  advance  by  original 
data  from  1500  hours,  8  October  1961.  For  comparison  Fig.  5  shows 
the  field  forecast  on  a  quasi-geostrophic  model,  and  in  Fig.  6 

the  actual  field  in  the  same  period. 

The  amount  of  mean  relative  error  was  calculated  for  19  points 
evenly  distributed  over  European  territory  of  the  USSR:  for  the 
quasi-geostrophic  forecast  for  the  given  period  e  *  0.75,  and  for 
the  quasl-solenoldal  -  t  •  0.68. 

Comparison  of  quasi-geostrophic  and  quasl-solenoldal  forecasts 
calculated  on  a  barotropic  model  for  the  northern  hemisphere  allows 
hoping  that  during  tests  of  quasl-solenoldal  forecasts  under 
operative  conditions  satisfactory  results  will  be  achieved. 
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Pig.  4.  H^qq 

cast  map  for  1500 
hours,  10  October 
1961  based  on 
original  data  for 
1500  hours, 

8  October  1961, 
calculated  accord¬ 
ing  to  a  quasi- 
solenoidal  model. 

A  continuous  line 
shows  the  contours 
of  the  area  for 
which  the  fore¬ 
cast  was  made. 


Pig.  5.  s500 

forecast  map  for 
1500  hours , 

10  October  1961 
according  to 
original  data  from 
1500  hours, 

8  October  1961, 
calculated  accord¬ 
ing  to  a  quasi- 
geostrophic  model. 


Fig.  6.  Actual  ^500  maP  **or  ^500  hours, 
10  October  1961. 
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THE  INVESTIGATION  OF  HEATING  IN  THE  STRATOSPHERE 
USING  NUMERICAL  EXPERIMENTS1 


B.  Barg  and  S.  A.  Uashkovich 


The  evolution  of  temperature  perturbation 
in  the  stratosphere  is  studied.  The  initial 
equations  were  proposed  by  A.  S.  Dubov  [1]  for 
forecasting  in  the  stratosphere.  The  qeustion 
abou1-  the  original  temperature  perturbation  is 
considered.  A  finite-difference  approximation 
of  equations  is  :.tated,  a  numerical  method  of 
solving  the  problem  is  formulated,  and  its 
calculating  stability  is  investigated.  Results 
of  calculations  for  different  variants  of  the 
initial  temperature  perturbation  are  presented. 


From  liter*  e  it  is  known  that  at  the  end  of  winter  as  well 
as  In  the  beginning  of  sprint  in  the  stratosphere  there  can  be  sudden 
warming  trends.  This  phenomenon,  discovered  by  Sgherhag  [10]  and 
known  by  the  name  "Berlin  phenomenon,"  was  then  confirmed  by  further 
observations  and  became  the  object  of  numerous  Investigations  (see, 
for  example,  [6,  9J). 


The  majority  of  investigations  were  in  first  place  to  establish 
possible  reasons  for  these  warming  trends.  This  work  does  not 
consider  these  questions  and  primary  attention  is  given  an  attempt 


‘The  present  article  is  the  result  of  the  joint  work  of  the 
Institute  for  the  Study  of  Large-Scale  Weather  of  the  GDR  Weather 
Service  and  Hydrometeorologic  Center  of  the  USSR. 
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to  study  the  evolution  of  the  warming  trend  and  its  values..  In 
other  words,  it  assumes  that  temperature  perturbation  in  the 
stratosphere  is  suitably  described  by  a  prlore  methods,  and  its 
further  development  requires  the  use  of  hydrodynamic  equations. 
Outwardly  such  an  investigation  is  similar  with  the  work  of 
Khlnkel'man  [4],  Berkovskiy  and  Shapiro  [73. 

However,  in  contrast  with  the  investigations  in  [73  we  limit 
ourselves  to  an  examination  of  processes  in  the  stratosphere. 
Therefore  we  can  use  the  work  of  Dubov  [13,  in  which  approximate 
forecast  equations  for  layers  of  the  atmosphere  with  a  high 
statistical  stability  are  obtained,  l.e.,  equations  which  are 
suitable  for  dynamic  examination  of  stratospheric  processes.  One 
ought  likewise  to  show  that  results  obtained  in  [13  allow  suffi¬ 
ciently  accurate  passing  from  equations  in  a  three-dimensional  space 
to  equations  for  a  two-dimensibnal  case.  This  circumstance  from  a 
mathematical  point  of  view  is  essential,  since  it  noticeably  facili¬ 
tates  calculations. 

According  to  Dubov  [13,  we  will  begin  from  the  following 
equations : 

■Hr— f<*  (i) 

Hr-.-fJm  4*>+<‘-  **>+-£.-£-]•  (2) 

Here  and  subsequently  we  accept  the  designations :  z  -  height  of 
considered  isobarlc  surface,  T  -  temperature  on  this  surface,  g  - 
acceleration  of  gravity,  l  -  Coriolis  parameter,  8  -  Rossby 
parameter,  A  -  Laplace  operator,  (c,  d)  -  Jacoby  operator,  v  -  dry 

d. 

adiabatic  gradient,  a  -  radius  of  earth. 

Initial  Conditions 

If  equations  (1),  (2)  are  used  to  solve  the  model  problem,  then 
it  Is  hardly  expedient  to  use  as  initial  conditions  the  data  of 
observations  for  practical  synoptic  situations.  It  is  more  useful 
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to  use  a  problem  which  gives  the  initial  fields  of  meteorological 
elements  in  the  simplest  form.  How  must  the  basic  fields,  and 
consequently,  the  initial  conditions  look? 

A  good  correspondence  to  actual  conditions  would  hardly  be 
achieved  if  mean  maps,  for  example  mean  monthly  maps,  would  be  used 
as  these  fields.  It  is  known  that  in  averaging  Important  parts  of  a 
field  are  lost,  such  as  the  intensity  of  stratospheric  polar  eddy. 

It  is  possible  to  approximate  more  the  practical  relationships  in  the 
atmosphere  if  we  suitably  stylize  the  determined  observable  typical 
states. 

To  do  this  the  chart  shows  profiles  of  the  geopotential  and 
temperatures  for  50  mbar  along  a  definite  meridian  at  0000  hours 
Greenwich  time  for  5,  16  and  17  January  1958  (Fig.  1). 


pig  1.  Approximation  of  H^q  meridian  profile  with  the 
aid  of  formula  (6).  1  -  meridian  distribution  by  data  of 

observations;  2  -  distribution  calculated  by  formula  (6). 
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Let  us  present  now  the  observed  profiles  using  an  analytic 
expression.  If  we  do  not  require  that  the  approximation  be  the 
best,  then  there  are  definite  degrees  of  freedom  which  one  could 
use. 

Prom  observations  which  must  generally  describe  the  polar 
vortex,  it  is  possible  to  borrow  the  following  facts: 

a.  At  the  pole  wind  velocity  is  zero. 

b.  Wind  velocity  turns  into  zero  between  the  equator  and  30° 

\ 

north  latitude. 

c.  Maximum  wind  velocity  is  between  60  and  75°  north  latitude 
(in  our  formulas  it  will  be  accepted  that  the  maximum  lies  at  60° 
north  latitude). 

/ 

d.  Undisturbed  wind  velocity  has  only  zonal  component  u. 

These  assumptions,  and  also  the  assumption  that  the  initial 
undisturbed  field  depends  on  longitude  X,  allow  selecting  for  wind 
velocity  a  parabolic  profile 

u~c  (»-*,)(»-»,),  (3) 

where  u  is  the  zonal  wind  velocity,  ft  -  polar  angle, 

8,  =  60C«2W;  <Uc=30e. 


From  expression  (3)  follows 

araC)»»3-C,9;  C)  = - 

*>MKC  *M»« 

Applying  the  formula  for  geostrophic  wind 


u  «= 


1  dz 
a  d  8  ’ 


(4) 


(5) 


we  obtain  from  ( 4 )  the  expression  for  the  height  of  the  50  mbar 
isobaric  surface 
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( z  is  expressed  in  decameters,  2q  ■  820  dam). 

In  Pig.  1  it  is  evident  that  the  selected  parabolical  meridian 
wind  profile  will  give  a  distribution  of  H j-q  which  agrees  well  with 
the  data  of  observations.  Such  means  can  also  obtain  the  expression 
for  approximation  of  temperature.  From  examination  of  the  field  it 
is  evident  that  the  investigation  must  embrace  the  range  from  the 
pole  to  rt=6ft#.  Accordingly  the  temperature  field  we  will  approximate 
also  in  the  zone  between  ti=0~  and  0"=60;;  in  this  ease  the  formula 
for  temperature  distribution  has  the  form 

T  - to  +  A  cos  3,56  +  273°.  (7) 

where  =  -70°,  A  =  12°.  As  can  be  seen  from  Pig.  2,  this  approxi¬ 
mation  is  sufficiently  accurate. 


These  expressions  approximate  undisturbed  basic  fields.  In 
accordance  with  our  original  thesis  we  assume  that  local  temperature 
perturbation  can  be  arbitrary.  Prom  systematical  considerations  before 
we  describe  how  this  perturbation  is  assigned,  we  will  give  certain 
Information  from  the  technology  of  a  numerical  solution.  On  a  map 
of  the  northern  hemisphere  in  stereographic  projection  (with  the 
section  covering  60°  north  latitude)  is  imposed  a  regular  grid  with 
42  x  42  points  (Fig.  3).  The  pole  does  not  coincide  with  the  grid 
point,  but  is  in  the  central  point  of  the  grid.  Grid  interval  s 
is  330  km  at  60°  north  latitude.  This  regular  grid  Is  used  as  the 
system  of  coordinates.  In  this  instance  temperature  perturbation 
T*  can  be  given  by  the  formula 


T*=A*  cos  + 


Here  x  ,  y  -  coordinates  of  the  center  of  the  temperature 
8  8 

perturbation,  expressed  in  units  of  length  equal  to  the  grid  interval; 
4*,  B*j  K*  -  arbitrary  constants,  r*  -  radius  of  temperature 
perturbation,  likewise  expressed  in  grid  intervals. 
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Prom  formula  (8)  It  Is  evident  that  perturbation  has  a  circular 
formula  relative  to  a  regular  grid  (on  the  spherical  earth  it 
acquired  a  pear  shape).  The  perturbation  lies  so  that  it  would  be 
as  far  as  possible  from  the  edges  of  the  considered  area  on  which 
artificial  boundary  conditions  are  assigned  (3s/9t  or  3T/9t  equal 
to  zero). 

Approximation  of  Derivatives  in  Time. 

Stability  Analysis  of  Solution 

When  selecting  a  numerical  method  for  solving  the  problem  at 
hand,  an  important  question  is  the  finite-difference  approximation 
of  derivatives  in  time.  Therefore  it  is  useful  to  estimate  the 
accuracy  of  the  solution  with  different  means  of  approximation  of 
these  derivatives.  Such  estimates  can  be  obtained  by  comparing 
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an  accurate  solution  with  the  numerical.  For  the  above  problem  it 
is  possible  to  find  an  accurate  solution  in  the  case  of  initial 
conditions  of  a  special  form.  Below  is  stated  the  means  of  obtaining 
the  accurate  solution  and  this  accurate  solution  is  compared  with 
the  numerical  with  different  variants  of  approximation  of  the 
derivative  in  time. 

We  will  consider  quasi-solenoidal  motion  and  assume  that  stream 
function  ^  is  connected  to  geopotential  a  by  the  following  approxi¬ 
mate  relationship: 


(9) 


Then  the  system  of  equations  of  our  problem  in  spherical 
coordinates  is  written  in  the  form 


(t,  n. 


(10) 


where  T  »  IT*. 

We  assume  that  the  stream  function  and  temperature  function 
can  be  presented  as  the  sum  of  two  components:  zonal  (not  depending 
upon  geographical  longitude  and  time)  and  nonzonal.  Then  we  write 
that 


tl*.  K  s  *)-*(•.  K  r*.  t); 

(11) 

Zonal  component  ijJ  we  represent  by  the  formula 

?— wg-^ga’co**.  (12) 

whe~e  a  -  index  of  circulation.  Then  from  the  equation  of  statics 
and  of  relationships  (9)  we  obtain  the  following  expression  for  the 
zonal  component  of  temperature: 


^7 


(13) 


Writing  T*  for  the  level  £?  *  C-^)  of  interest,  we  obtain 

7H-1.—  Cl- *»co«*,  (1*0 

where 

c •  “  -^(4$%/  **• 

Let  us  substitute  now  relationships  (11),  (12),  (14)  In  equations 
(10): 


d*T*  c,rn  >  tv.  ,  .W 

+  A*'>- 


-k  *T*  ■ 

r 


Ai  (u  —  tWi  »*ii  * 


(15) 


where 


*,-*«  +  .)  + 

(t.  - 

k  _  (bYt  ■ 
** 


We  accept  subsequently  that  and  can  be  approximately 

considered  constants.  ' 

We  assume  that  in  initial  moment  t  s  0  the  nonsonal  parts  of 
stream  function  ip'  and  temperature  T*’  are  represented  by  the 
formulas : 


f  U-Re  VW). 
7-*'U-Re*e*,"*/3r<»). 


(16) 


where  P ™  is  an  adjoint  Legendre  polynomial. 


Accordingly  look  for  the  solution  of  system  (15)  in  the  form: 
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Lee  us  go  now  to  the  solution  of  system  (18)  by  the  numerical 
method.  Namely,  let  us  examine  two  variants  of  approximation: 

1)  explicit  scheme  with  replacement  of  derivatives  in  time  by 
"forward”  differences,  2)  implicit  scheme. 

In  the  first  variant  derivatives  in  time  are  represented  in 
the  form: 


4A  /('“-At  48  -a* 

4t  “  it  41  “  it 


(23) 


and  all  remaining  terms  in  equations  ( 1 8 )  are  written  for  preceding 
instant  t.  We  set 


A -AS'1;  B-BS' 


(2H) 


Then  considering  expressions  (23)  and  (2^4)  equations  ( 1 8 )  will 
take  the  form: 


** •*'  <(  -$(2  -  y)A,  +  [>  y  -  2(.  +  -)!*,). 


(25) 


From  (25)  can  be  obtained  the  following  equation  for 
determination  of  B^/K^\ 


(_*  j_ +  -£<2 _ y) _ o. 


(26) 


The  quantity  B^ /A^  can  be  both  real  and  complex.  Let  us  write 
down  the  expression  for  B 2^2  *n  t^ie  ^ortn 


&)-%  =  /*  o. 


(27) 


while 


/ 1  when  /,*0. 


Substituting  (27)  in  (25)  and  separating  real  and  imaginary 
iarts,  we  obtain  after  simple  transformations: 
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y  ±m\ 


(28) 

(29) 

(30) 


It  is  obvious  that  when 


01 

we  have  v,<o  and,  consequently,  in  accordance  with  (29),  (30)  ampli¬ 
tude  grows  in  time. 


Let  us  examine  now  the  solution  in  the  case  of  the  implicit 
scheme.  In  this  variant  derivatives  in  time  also  are  approximated 
by  the  formula  (23),  however  all  remaining  terms  of  (18)  are  written 
for  the  subsequent  moment  of  time  t  *  fit.  Using  (23)  and  (24),  we 
can  write  system  (18)  in  the  form: 

fy  -  lm(k,  -  *  y)8  t]Aj  +  lm*t  3  tB,  - 

\y-\-lmlt\ay  —  2(i  -f  *)|  Jfl,  -  ^  (2  -  y)lm  8  tA,  -  yB*-^'  (  31 ) 

Prom  (31)  we  obtain  the  equation  for  the  determination  of 
B2/A2t  completely  coinciding  with  (26).  Making  use  of  expression 
(27),  after  simple  transformations 


Obtained  formulas  allow  evaluating  the  accuracy  of  the  numerical 
solution.  For  this  let  us  compare  results  of  calculation  for 
accurate  solution  of  (17)  and  (19)  with  calculations  from  formulas 
(28)  and  (32). 


Calculations  were  produoed  for  the  following  values  of 
parameters:  fit  ■  1  hour,  3*^  ■  200° ,  -  250°,  y  ■  0,  m  «  1,  «  ■  2, 

5,  10,  Very  important  is  the  seleotion  of  the  values  of  circulation 
index  and  its  derivative  in  the  vertioal  coordinate.  The  last 
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quantity  is  substantially  different  in  different  seasons.  Inasmuch 
as  sudden  warming  trends  in  the  stratosphere  are  observed  predomi¬ 
nantly  at  the  end  of  winter  and  beginning  of  spring,  the  value 
dcx/dz,  referred  to  this  time  of  year  was  taken,  namely 

Results  of  calculations  are  presented  in  Table  1.  The  table 
shows  quantities  characterizing  the  change  in  wave  amplitude  in  the 

time  interval  6t  =  1  hour.  In  the  case  of  an  accurate  solution  there 

(/) 

exist  only  real  roots  of  characteristic  equation  (21),  i.e., 
therefore  wave  amplitude  does  not  change  in  time.  With  the  numerical 
solution  method  there  appear  changes  of  amplitude  in  time.  These 
changes  seem  greatest  for  small  values  of  n,  so  for  n  =  2  in  the 
hourly  interval  amplitude  changes  by  0.3^  of  the  initial  amount. 

Using  "forward"  difference^  there  is  a  growth  of  amplitude;  this 
variant  of  solution  seems  unstable  in  calculations.  In  the  case  of 
the  implicit  scheme,  stable  in  calculations,  a  decrease  of  amplitude 
in  time  is  noted. 


Table  1 . 


Ace urate 

lolution 

J  "forward"  difference 

Implioit  eoheme 

n 

(1) 

: 

exp 

exp 

i  .  . 

2 

1.0 

1.0 

1.0030 

1,0014 

0.9970 

0,9986" 

5 

1.0- 

1.0 

1,0006 

1.0000 

0,9992 

1,0000 

10 

.  1.0 

ID 

1.0)01 

1.000 

•  0,999$  " 

1,0000 

Calculation  System 


Passage  to  a  Cartesian  system  of  coordinates  is  carried  out. 
A  scale  factor  is  introduced  by  the  formula 


_  «<!  +  ►»» 
a  • 


(33) 
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where  c-a(l+sin$o):p»»-~;  f  -  geographical  latitude  on  which  the  plane 
of  projection  lies;  r  -  distance  of  considered  point  from  the  pole  in 
the  plane  of  projection. 

The  Laplace  operator  in  the  right  sides  of  (1)  and  (2)  is 
approximated  by  the  expression 

25t[*1  +  +  Z*  +  *4  +  -J-  (*i  +  +  *»)  —  6*0  J.  ( 3  4 ) 

The  arrangement  of  points  is  shown  in  Pig.  4.  Terms  in  the 

Jacobian  are  represented  in  the  form 

•af  if  *-sr[“«  -  “» + 2<a|  -  *1) + **■- *•]  L  + -  •  *Jc  (35) 

We  introduce  the  designations: 

(*.  p).  (36) 

The  expression  of  the  Laplace  operator  in  the  left  sides  of  (1) 

and  (2)  we  approximate  by  the  formulas: 

**  7T  xr^£i + *  ** + — **  -7^7 

(37) 

while 

di  ^  i*  dT  ^  &  r 
~5T~~  it  ’  Tt  it  • 

We  will  use  the  MTS  system  for  all  quantities  except  geopotential  z; 
we  will  express  it  in  decameters.  Designate  l  =  lQ\)  then  equations 
(1)  and  (2)  can  be  represented  in  the  form: 


i  •  -> 

to 


(38) 


(39) 


+  M  '?  *!»  (*••?• 

86400 

where  ,  n  -  number  of  time  intervals  in  a  twenty-four  hour 

period. 
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Fig.  4.  Arrangement 
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The  solution  of  equations  proceeds  as  follows: 

a.  The  right  sides  of  equation  (38)  are  resolved. 

b.  Determine  6z  in  decameters  by  solving  the  Poisson  equations 
and  calculate  z  for  tQ  +  6t.  As  boundary  conditions  with  the  solution 
of  the  Poisson  equation  we  take  te|«p-*0  (first  boundary  problem  or 
Dirlchlet  problem). 

c.  Calculate  the  right  side  of  equation  (39)  using 

d.  By  means  of  the  colutions  of  the  Poisson  equation  with 
boundary  condition  *7* 0  determine  6T  and  find  T |»,+w  (then 
everything  is  repeated). 

Solution  of  Poisson  Equations 

It  is  known  that  in  solving  forecasting  problems  which  are  close 
to  ours  the  use  of  the  Poisson  equation  gives  results  which  are 
physically  insufficiently  satisfactory.  In  order  to  understand  the 
difficulties  it  is  necessary  to  take  into  account  that  in  the 
considered  case  it  is  necessary  to  solve  an  elliptical  equation 
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with  a  known  right  side.  In  this  instance  the  effect  of  perturbation 
is  extended  with  an  infinitely  high  rate  and  each  point  of  the  field 
instantly  feels  the  effect.  Prom  the  viewpoint  of  meteorology  this 
leads  to  unrealistic  results.  In  order  to  avoid  these  effects,  the 
Poisson  equation  is  replaced  by  the  Helmholtz  equation.  The  basis 
of  such  a  replacement  is  different  arguments,  which  bear  partially  a 
physical  character,  partially  a  mathematical.  A  survey  of  these 
questions  can  be  found  in  [5]. 

In  a  more  of  less  arbitrary  passage  to  the  Helmholtz  equation, 
written  in  the  form 

lA-WIt-F  (HO) 

(A  in  dimensionless  variables),  the  question  arises  about  determina- 

p 

tion  of  (sfe)  .  There  are  no  definite  rules  on  this  and  the  question 
can  be  solved  only  experimentally.  However,  it  is  not  clear  what  the 
results  obtained  using  the  solution  of  the  Helmholtz  equation  (HO) 
with  an  assigned  value  of  (efe)  should  be  compared  with. 

To  solve  equation  (HO)  we  use  a  method  proposed  by  S.  V. 
Nemchinov  [3]  for  the  case  of  a  right-angled  grid.  The  advantages 
of  this  method  are  obvious:  on  one  hand,  because  of  the  use  of 
recurrent  formulas  a  minimum  number  of  memory  cells  is  necessary, 
on  the  other  a  relatively  small  number  requires  arithmetical 
operations.  To  this  one  ought  to  add  that  the  solution  has  a  high 
degree  of  accuracy,  which  when  using  iteration  methods  can  be 
achieved  only  with  great  difficulty,  and  that  the  calculating  process 
is  stable  with  respect  to  rounding  error.  The  solution  with  a  change 
in  fe  does  not  require  a  change  in  the  program  (in  the  case  of  k  =  0 
the  solution  cf  the  Poisson  equation  is  automatically  obtained). 

In  the  practice  numerical  forecasts  are  satisfied  by  attempts 

n 

to  obtain  an  optimum  forecast  and  in  this  manner  determine  (afe)  . 

In  order  to  use  this  experiment  in  this  case  a  special  Investigation 
would  be  necessary.  Therefore  we  used  another  way. 
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Several  years  ago  S.  L.  Belousov  [2]  proposed  a  means  of 
solving  the  Poisson  equation,  which  since  then  has  been  successfully 
used.  The  formulas  are  best  of  all  called  a  local  solution;  we  will 
use  this  term  subsequently. 

The  solution  of  the  Poisson  equation 

^ CD 

at.  known  boundary  values  can  be  presented  in  the  form 

fc— -if  (42) 

o  o  > 

where  -  desired  magnitude  of  the  function  in  the  center  of  a 
circle  of  radius  R. 

For  calculations  according  to  formula  (42)  Belousov  used  the 
method  of  series  approximations.  As  a  first  approximation  he  assumes 
that  the  contour  integral  in  (42)  is  absent,  and  to  calculate  the 
first  integral  in  (42)  he  used  the  formula1 

-  s2\~Fn  +  ±(F,  +  Ft  +  FtS-\-±[F,  +  ...  +f„)J.  (43) 

Since  the  quantity  i|/q  is  found  in  all  pchlnts,  the  approximate  value 
of  the  contour  integral  in  (42)  is  calculated.  In  the  contour 
integral  the  value  of  i|>q  is  substituted;  the  second  approximation 
Is  calculated  by  the  formula 

+  (44) 

i-l 

'if  we  consider  a  relatively  small  neighborhood  of  point  0,  then 
the  approximate  value  of  the  integral  can  be  obtained  by  approxi¬ 
mating  the  Laplace  operator  by  the  formula 

*  -r(*.  +.♦»  +  *4  -  <*>) . 

In  this  case  the  Laplace  operator  was  written  for  9  points  (from 
zero  to  8th),  and  the  values  of  i|>  at  points  from  the  9th  to  the  20th 
are  taken  as  zero.  The  solution  of  the  obtained  system  of  algebraic 
equations  leads  to  formula  (43). 
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We  used  this  means  of  solving  the  Poisson  equation.  Not  stopping 
here  on  questions  connected  with  the  use  of  this,  we  refer  the 
reader  to  [5,  8]. 

It  is  possible  to  compare  results  for  two  variants  of  the 
solution  of  system  (38),  (39).*  1)  using  the  Helmholtz  equation  for 

various  values  of  (sfe)  ;  2)  using  the  local  solution  considering 
two  approximations.  It  is  established  that  results  will  well  agree 
if  in  the  first  case  we  take  (ek)  =  0.4.  This  value  already  has 
been  widely  used  in  actual  numerical  forecasts.  Thus  the  local 
solution  can  be  successfully  used. 

Preliminary  Results 

Solving  the  problem  by  computer,  initially  a  finite-difference 
approximation  of  derivatives  in  time  with  the  aid  of  "forward" 
differences  was  used.  However',  the  calculations  showed  that  despite 
a  decrease  in  the  time  interval  t^,  30  minutes,  the  system  is  unstable. 
At  t  -  +  48  hours  stray  waves  show  up  and  intensify,  thus  hampering 

the  interpretation  of  results.  This  agrees  with  theoretical  estimates 
obtained  above.  Therefore  in  basic  variants  of  the  calculations  an 
implicit  scheme  was  used.  The  time  interval  in  this  case  could  be 
taken  as  3  hours. 

Calculations  were  carried  out  for  the  following  values  of 
parameters,  determining  temperature  perturbation:  A*  =  K*  =  6°; 

B  =  3*6;  r*  <  5.  In  this  case  cases  were  considered  when 
x  =  y  =  -7  and  x  -  y  =  -11.9.  The  basic  field  of  flows  was 
given  by  the  values  of  parameters: 

1)  C1  “  —  7T ;  ci "  x 30  m/s  ; 

2) c,  “  ~«r  ;  c*  “  tr;  *«•»« — 10  m/s  • 

On  the  basis  of  results  of  calculation  the  following  preliminary 
conclusions  can  be  made. 
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1.  A  temperature  perturbation  of  a  shown  intensity  can 
substantially  influence  circulation  in  the  stratosphere.  It  does  not 
lead  therefore  to  a  breakdown  of  the  polar  vortex. 

2.  As  a  result  of  temperature  perturbation  a  perturbation  in 
the  geopotential  field  is  developed,  the  amplitude  of  which  grows 
during  the  considered  time  interval,  while  the  amplitude  of  the 
temperature  perturbation  diminishes. 

3.  The  maximum  of  the  geopotential  perturbation  remains  during 
the  first  48  hours  at  the  same  latitude  on  which  was  the  temperature 
perturbation.  Then  in  all  cases  the  geopotential  perturbation  shifts 
to  the  north,  which  agrees  with  the  data  of  observations.  However, 
the  temperature  perturbation  remains  at  the  original  latitude.  The 
last  result  strongly  differs  from  processes  observed  in  nature.  This 
difference  is  partially  conditioned  by  simplifications  of  the  model, 
and  partially  by  neglecting  effects  connected  with  the  ozone. 

4.  The  shift  of  the  geopotential  perturbation  to  the  north  is 
expressed  stronger  the  better  the  polar  vortex  is  developed,  i.e., 
the  greater  the  meridional  wind  shift. 
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ON  THE  APPLICATION  OF  THE  FILTRATION  THEORY  OF  RANDOM 
PROCESSES  TO  SOME  PROBLEMS  OF  OBJECTIVE  ANALYSIS 


B.  V.  Ovchinskiy 

Formulas  for  calculating  weighting  functions 
of  optimum  interpolation  and  smoothing  meteorolog¬ 
ical  elements  if  correlation  functions  of  the 
"true"  quantity  (signal)  and  errors  of  observa¬ 
tions  (noise)  are  known. 


I .  Formulating  the  Problem.  Derivation 
of  Basic  Equations 

As  is  known  the  data  of  observations  contain  random  errors. 

During  optimum  Interpolation  it  is  necessary  to  consider  this  and 

/ 

treat  the  meteorological  element  at  an  interpolated  point  free  from 
these  random  errors.  '■ 

Even  when  the  field  of  meteorological  elements  is  presented 
rather  fully,  for  further  analysis  It  is  useful  to  remove  from  the 
field  unsystematic  errors,  i.e.,  smooth  the  field. 

As  has  been  accepted  in  objective  analysis,  a  field  of  meteoro¬ 
logical  elements  is  given  at  discrete  points  3^,  h z,  ...,  Each 

of  the  quantities  5^,  except  the  true  value  ,  still  has  errors  e_. , 
so  that 


(1) 
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Concerning  errors  e.,  they  can  In  the  simplest  case  be  either 
purely  random,  or  contain  small-scale  fluctuations,  of  the  "noise" 
type,  as  happens  In  radiophysics . 

In  objective  analysis  errors  were  taken  as  purely  random  and 
satisfied  the  following  conditions  [2]: 

1)  mathematical  expectation  (mean)  error  e.  is  equal  to  zero, 
i.e. ,  Eft  .)  •  0; 

X 

2)  errors  e.  are  not  correlated  with  true  value  of  the  meteoro- 

X 

logical  element  Eft. g.)  •  0; 

X  X 

3)  errors  t^  do  not  correlate  with  each  other,  i.e.,  Eft^tJ  = 

■  0  when  i  ft  j . 

If  we  treat  more  broadly  the  quantities  e .  and  include  here 
small-scale  fluctuations,  then  we  must  reject  the  assumption  about 
the  noncorrelatability  of  the  e..  In  meteorology  this  circumstance 
was  studied  by  Thompson  [8],  using  the  apparatus  of  the  theory  of 
filtration  of  random  processes  for  optimum  smoothing  of  meteoro¬ 
logical  fields. 

Thus  we  assume  that  g.  and  e.  are  stationary  random  quantities, 

X  X 

which  obey  conditions  1)  and  2).  The  correlation  functions  of  g. 
and  e  .  will  be : 

t 

where  t..  —  distance  between  i-th  and  j-th  stations. 

^  J 

In  relationship  to  correlation  functions  R  fr.J  and  R  fr.J 

g  e  tj 

let  us  note  that  radius  correlation  Is  considerably  less  than  the 
radius  of  correlation  of  the  g^. .  Tnis  mer.ns  that  the  connection 
between  the  e.  weakens  much  faster  with  distance  (or  with  tine), 
than  the  connection  between  the  g.. 
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In  accordance  with  the  theory  of  optimum  interpolation  we  put 
together  the  expression 

***£**-%*  +  '&  (2) 

and  select  weights  such  that  dispersion 

(3) 

is  minimum.  After  this  the  interpolated  value  of  the  meteorological 
element  will  not  contain  "noise,"  and  thus  is  a  smoothed  quantity. 

Let  us  suppose,  as  usual,  that  h .  are  standard  deviations  and, 

u 

according  to  the  first  condition  relative  to  errors  e.,  we  have 

0. 

After  this  formula  (3)  can  be  brought  to  the  form: 


P  -  £(*’)  -2|/?,(T0()P(+v|^tiy)  +|| R.h)P,P, .  C  *0 

J 

2 

We  divide  both  sides  of  equality  (4)  by  a  and  turn  to  standard 
ized  correlation  functions  y  .  .  for  the  signal  ("true  quantity")  and 

ij 

v..  for  errors  e.  ("noise"),  then  formula  (4)  can  be  rewritten 


HtP' + f . 1  + T  %'uP,Pj 


(4') 


2  2 

The  quantity  o/o  in  radiophysics  Is  called  the  signal-to- 

E  O 

noise  ratio.  Subsequently  we  will  designate  it  by  A.  For  minimum 

2  2 

J !  it  is  necessary  that  3c7^/3P^  =  0.  Hence  we  obtain  a  system  of 
linear  equations  for  determination  of  weights 

/- 1,2, . AT.  (5) 
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Assuming  that  errors  e.  do  not  correlate  among  themselves,  i.e., 

they  are  purely  random,  then  function  v.  .  can  be  represented  as: 

t  J 

f  I,  if  /-/; 

'-{lit  i+j.  (6) 

This  case  corresponds  to  the  fact  that  in  radiophysics  it  i3 
accepted  to  say  "white  noise"  [6],  For  white  noise  equation  (5) 
can  be  rewritten  in  the  form: 


..  (5') 

Formula  (5*)  coincides  with  the  usual  equation  used  in  objective 

2  2 

analysis  [2],  because  n  =  X  =  a  /a  .  A  work  of  Thompson  [9] 

£  8 

expresses  certain  considerations  about  X.  The  signal-to-noise  ratio 
fluctuates  from  X  =  0.1  to  =  0.33. 

If  weights  P.  are  determined  according  to  formula  (5),  then 

‘V 

interpolation  error,  according  to  (4*),  can  be  found  by  the  formula 

(7) 

i-i 

Let  us  note  that  formula  (7)  contains  the  correlation  function 
of  noise  v^..  In  order  to  explain  in  general  terms  the  effect  of 
the  noise  part  of  the  field  on  interpolation  error,  we  set  X  =  1 
and  y..  =  v...  Then  from  equation  (5)  we  obtain 

tj 

i.e.,  weights  P.  will  be  half  as  large  as  in  the  same  "noiseless" 

J 

field,  and  from  formula  (7)  it  follows  that  interpolation  error 

.  .v 

increases  and  will  be  **=d.v  PiHi-  This  increase  in  interpolation 

error  for  such  a  noise  field  can  be  explained  by  the  fact  that  there 
are  two  identical  superimposed  fields  here. 
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Let  us  turn  again  to  system  of  equations  (5)  for  determination 
of  weights  of  optimum  interpolation.  Let  us  write  this  system  for 
the  case  wnen  the  data  of  observations  have  been  taken  through  equal 
intervals  in  3pace.1  This  simplified  scheme  will  help  us  to  investi¬ 
gate  in  more  detail  the  structure  and  features  of  the  weighting 
function,  which  is  the  basic  task  of  this  work.  Below  we  will 
examine  the  case  when  observations  of  the  quantities  in  meteorolog¬ 
ical  elements  are  continuous.  Thus  let  us  assume  that  the  quantities 
in  meteorological  elements  3^  are  given  through  identical 

distances  from  one  another,  and  for  simplicity  let  us  take  these 
distances  as  one.  Then  the  correlation  function  can  be  written: 

£<&*/) 

For  standardized  correlation  functions  let  us  introduce  new 
designations : 

After  this  system  (5)  for  determination  of  weights  of  optimum 
interpolation  can  be  replaced  by  other  system  of  equations 

“>)+*■  rAk-j)]P,-rJik\  k»  lt  2,  t . N.  (8) 

Equation  (8)  in  statistical  dynamics  of  pulse  systems  [5]  is 
called  a  discrete  analog  the  Wiener-Hcpf  equation. 

Let  us  assume  that  assa  result  of  observations  all  possible 
values  of  3^  have  been  obtained  including  Our  problem  as  before 

is  the  determination  of  gQ  on  the  basis  of  the  observations.  In 

*0f  course,  we  need  not  hold  strictly  the  selection  of  observa¬ 
tion  points  through  equal  intervals.  Deviations  are  possible  in 
sucn  quantities  which  will  not  substantially  affect  the  value  o** 
the  correlation  function. 
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this  instance  we  say  that  we  smooth  It Q .  Formula  (2)  will  have  the 
form: 


-  gw  +«|) . 


(9) 


The  distinction  from  formula  (2)  is  that  the  summing  begins 
from  i  »  0. 

System  of  equations  (5)  for  determination  of  weights  P.  remains 
as  before  with  the  difference  that  the  summing  must  begin  with  i  =  0; 
3=0.  M.  I.  Yudin  proposed  two  formulas  for  smoothing  meteorolog¬ 
ical  fields,  which  in  form  correspond  to  expression  (9)  and  are 
written  thus  [10]: 

,  :s  ,r :C 

H  .  0,36/*,  +  0,06  2  +  0,06 (//,  +  Hn) . 

-  ^ 
r 

The  stations  with  the  numbers  1-6  lay  along  the  circumference  with 
center  and  stations  8,  11  were  outside  the  circumference  symmet¬ 
rical  to  the  origin  of  coordinates. 

Let  us  examine  the  problem  of  smoothing  meteorological  elements 

% 

u(x) ,  when  x  takes  continuous  values.  Of  course,  for  discrete 
values  of  x  the  quantity  u(x)  must  coincide  with  the  observed  quan¬ 
tities  lt1,  1t2,  a^,  i.e.,  u[0]  »  3*0,  u[l]  =  w[y]  = 

-  It 


Just  as  earlier  we  assume  u(x) ,  except  the  true  value  (signal) 
g(x),  still  contains  random  error  (noise)  t(x) .  Furthermore, 

u(x)—g(x) +*(*)• 

The  smoothed  out  value  of  u(0)  must  not  contain  noise  and  its 
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approximate  magnitude  by  analogy  with  (8)  will  be1 

m 

go^ju(x)f(x)dx.  (!0) 

Let  us  find  such  weight  function  f(x)  that  the  mean  quadratic 
error  of  approximate  equality  (10)  is  minimum  [11] 

P(f)  -*(*.-  [  !*(■*)  +  *(x)]f(x)dx  }S  ( 11 ) 

We  transform  preliminarily  formula  (11),  making  use  of  the  fact 
that  g(x)  and  e(x)  are  stationary  random  functions,  not  correlated 
with  one  another,  i.e.,  £[e(xj  g(x)l  -  0.  As  a  result  we  obtain 

A/)-  ^0)  -  2 1-  Rt(x)f{x)dx  +jj  \Rg(x  -  y)  + 

+  rt,(*  -  y)\f(x)/(y)dxdy.  ( 12 ) 

o 

In  order  to  find  the  minimum  functional  J  (f) ,  it  is  necessary 
in  formula  (12)  instead  of  f(x)  to  substitute  f(x)  +  y$(x) .  The 
necessary  condition  of  extreme  of  the  functional  will  be  [3] 


Having  completed  the  calculations,  we  pass  to  the  integral  equation 
for  determination  of  f(x) 

j  -  y)  +  RM  -  y)\f(y)dy  -  Rjx).  (13) 

It  is  possible  to  show  a  sufficiency  of  condition  (13),  in 

0 

other  words,  if  /(*)  satifies  integral  equation  (13),  then  J  (f) 
will  be  minimum  [7].  Equation  (13)  is  widely  used  in  the  solution 
of  problems  of  automatic  control  and  is  called  the  Wiener-Hopf 

*The  upper  limit  of  the  integral  is  equal  to  infinity,  which  is 
not  obligatory.  However,  in  the  beginning  we  take  this  for  the 
sake  of  simplicity  of  further  formulas. 
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integral  equation.  Because  usually  it  is  more  convenient  to  deal 
with  normalized  correlation  functions,  equation  (13)  we  rewrite  in 
the  following  form: 

J  M* -y)-H-  r.  ix  -  y)J  f(y\dy  -  r£x) .  ( 1 4 ) 

Equation  (14)  is  the  original  for  determination  of  weight 
functions  of  smoothing. 


Weight  Functions  of  Optimum  Interpolation 
and  Smoothing 


Let  us  pass  now  to  the  solution  of  system  (8) 


2  W*  -/) + *  r.(k  -j)]Pj  -  rtfl  *-»  I,  2,  3,  ...  ,  JV. 


The  means  of  solving  equation  (8),  which  is  presented  later, 
was  used  in  statistical  dynamics  of  pulse  systems  [5]  as  well  as  in 
problems  of  automatic  control  [1,  7]. 


By  approximating  the  correlation  function  by  the  sum  of  model 
functions  R(t)  *  V  we  can  find  the  P.  solutions  by  the  method 

Ci  1 

of  indeterminate  coefficients.  Subsequently  we  will  limit  ourselves 
to  approximating  the  correlation  function  with  only  the  first  term 
of  the  sum,  i.e.,  we  set 

(15) 

where  a  is  determined  according  to  empirical  correlation  function. 


So  for  the  two  meteorological  elements  by  which  we  will  conduct 
the  calculation  (wind  and  dew  point),  the  distinctions  of  the 
approximated  curve  from  the  observed  correlation  coefficients 
(Figs.  1  and  2)  we  consider  as  entirely  permissible.  Therefore  we 
assume  that  the  correlation  function  for  the  dew  point  will  be 
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Pig.  1.  Autocorrelation  posi¬ 
tion  function  of  the  dew  point 
(spring)  [2].  1  —  observed, 

2  -  according  to  the  formula. 


Pig.  2.  Autocorrelation 
function  of  wind  [4].  1  - 

observed,  2  -  according  to 
the  formula. 


rg(p)  mm *-«•«* ,  and  for  wind1  r((p)  —  *-*•**.  where  distance  p  is  expressed 
in  arbitrary  units,  shown  in  Figs.  1  and  2.  Let  us  examine  from  the 
beginning  the  case  when  the  interference  is  white  noise.  The  corre¬ 
lation  function  for  errors  [see  formula  (6)]  can  be  written: 


during  calculations  it  was  accepted  r^i-#-****. 


I.  if  *-/; 
0,  if  k+j. 


Systems  of  equations  for  determination  of  weights  during  optimum 
interpolation  and  smoothing  can  be  represented  in  the  form: 


+  >P*-r<(*)t  *-l,  2,  /;/.  N.  (16) 

r/jt),  k-0,  I,  2,  N.  (17) 

Let  us  examine  from  the  beginning  system  (16).  Its  solution  we  look 
for  in  the  form1 

pj  =  Ae-',+Be'1.  «  (18) 


Coefficients  A,  B,  y  will  be  determined.  Let  us  substitute 
expressions  (15)  and  (18)  in  original  system  (16),  then  we  obtain 


2 ]iAe-'J+Be'1]  1 + Be' /)e-u"*)+ 

or  f 

Be-  e(T+,u  + 

+  fle‘*£  «<7'*w  +  >•  Ae~"+  /.  «T  * « e~* *,  A  - 1,  2,'...  ,  JV. 
>-*+ 1 


Having  completed  the  summation  and  equating  coefficients  In  the 
right  and  left  sides  at  e_Y^(eY^),  eY^  we  arrive  at  the  equa¬ 

tions  for  determination  of  y ,  A,  B: 


WT  -  #*  H#1  -  #-•) 


0. 


At'  Bt-' 

f’-i1  "’"ru/  "* 


I, 


A*''*  . 

»i-i-  t  r’-i- 


(19) 

(20) 


Selection  of  the  form  of  the  solution  will  be  clarified  at  the 
end  of  the  article. 


69 


The  order  of  solution  of  the  system  is:  from  equation  (19)  we 
determine  y»  an d  then  from  system  (20)  we  find  A  and  B.  For  the 
calculation  of  y,  A  and  B  it  is  possible  to  make  use  of  the  following 


formulas 


where 


f,  imi + 18— n  +  yV— i  T  mi +P8)  —  w/*' 

7kp 


.  Ll  1  )»(«'  -«*) 
ft 

+  (#T-  AV1 


From  formula  (21)  it  follows  that  y  depends  on  (quantity  of 
stations).  The  second  root  of  equation  (19)  gives  a  negative  value 
of  y,  which  does  not  have  physical  meaning.  Figure  3  shows  y  as  a 
function  of  X. 


A 


4?  «♦  5#  <0* 


£  5  a?  55  5 *c  *s sox 

“tg.  3.  y  as  a  function  of  X. 


As  it  appears  in  the  graph,  an  increase  of  X  decreases  values 
y  as  they  asymptotically  approach  a.  However,  this  asymptotic 
behavior  is  achieved  for  X  of  the  order  of  AO-50,  which  is  unreal¬ 
istic.  Therefore  it  is  possible  to  assume  that  always  y  is  more 
than  a.  In  other  words,  the  index  of  weight  function  y  is  much  more 
than  the  same  index  for  the  correlation  function. 


The  turn  now  to  determination  of  the  weight  functions  of  optimum 
smoothing.  For  this  we  must  solve  system  of  equations  (17). 

The  weight  functions  we  look  for  (as  earlier)  in  the  following 

form: 

The  solution  of  system  (17)  follows  exactly  system  (16). 

As  a  result  of  this  solution  we  obtain  y  «  y^ »  i.e.,  the  index 
of  damping  of  the  weight  function  for  optimum  y  interpolation  coin¬ 
cides  with  the  same  index  y^  for  smoothing  and  is  determined  by 
formula  (21).  Coefficients  A^  and  B 1  can  be  found  by  the  formula: 


.  _ (#»-•-!)»  _ 

«*)(!—  #-«7— «»|  * 


(2*0 

(25) 


Finally  let  us  examine  the  case  when  It  (number  of  stations) 
increases  without  limit  (it  *).  Then  from  formulas  (23)  and  (25) 
it  is  evident  that 


B-+  0  and  B ,  -+Q. 

The  weight  functions  for  optimum  P.  interpolation  and  smoothing 

J 

have  the  form: 

j- o.  1. 1, ... 

The  change  in  weight  function  depending  on  A  for  optimum  inter¬ 
polation  P.  *  A and  smoothing  £ .  «  A is  shown  in  Fig.  b . 

On  the  graphs  are  values  of  weights  only  for  P PQ.  By  these 

data  it  is  possible  to  establish  that  weight  ^  is  more  than  weight 
Pj  for  all  X. 
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Pig.  Change  in  weight 
function  depending  on  X  for 
smoothing  wind  (1)  and  dew 
points  (3)  for  optimum  interpo¬ 
lation  of  wind  (2)  and  dew 
points  (*0. 


Let  us  turn  now  to  the  more  general  case  when  errors  e.  are 
connected  between  one  another  and  the  correlation  function  is 
assigned  by  formula 

rjft-e-H*. 

Furthermore,  we  assume  the  upper  limit  of  summation  N  is  infi¬ 
nite;  then  system  (8)  can  be  written 


V  \rj*  —/)  +  >•  rXk  -  j)\P ,  -  rjkh  k  -  I.  2.  3 .  ( 8  -  ) 

i-i 


The  solution  of  system  (8')  we  look  for  in  the  form: 

P,-Ce «~0,\j-\].  /-I.  2.  ...  (26) 


Quantities  d,  C,  D  will  be  determined;  o[j]  is  the  pulse  single 
function  which  is  assigned 

Ul  I  0.  lf  J  +  0. 

Let  us  substitute  in  both  sides  of  equation  (8')  instead  of 

r  (x).  r  (x)  and  P.  their  expressions,  then  we  arrive  at 
g  E  3 
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v  [*-**-/> + \Ct-“  +  DU-UI  +  i 

f?i  f&t 

+).e-»>-"]Ce~* **•!.*.  ... 


Summing  up  and  equating  coefficients  in  the  right  and  left 
sides  at  ,  a~a^t  we  arrive  at  equations  for  determination 

of  d,  C  and  D : 


e'  —  e 


•  +  *•  ■ 


r  —  t 


C  _■ 

^  c 

e'  -  t>.  - 


+  o 

+  D 


0. 


-0. 


(27) 

(28) 


If  we  must  smooth  It. ,  then,  as  already  was  noted  above  in 

formula  (8'),  the  summing  must  begin  with  j  -  0.  The  weight  function 

we  find  according  to  formula 

\ 

r 

where  dj  as  before  is  determined  from  equation  (27);  and  can 
be  found  by  solving  system  of  equations: 


c, 

i 

C, 

i  - 


+  P, 


I. 

0. 


(29) 


Hence  it  follows  that  the  difference  in  weight  functions  of 
Interpolation  and  smoothing  is  only  because  of  coefficients  C  and  D. 


The  index  of  damping  of  weight  function  d  can  be  found  by 
solving  equation  (27) 


-  do  ihi  +  r)  , 

ij7>  -  w  ^  r-  + 

-  nr*  -mi  -  ^  )i- -  <Ke 

-  Dr*'  -  Up*  -  i)H 


np- 


1)^ 


(30) 
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It  is  accepted  that  0  -  m a  and,  as  noted  above,  0  >  a,  so  that  m  >  1. 

Figures  5  and  6  show  changes  of  e ^  depending  on  X  at  m  =  2, 
m  *  3  for  th£  dew  point  and  wind.  For  large  X  the  quantity  d  tends 
to  asymptotically  to  a.  The  difference  between  a  and  d  becomes 
unimportant  beginning  with  X  *»  25  and  more  for  different  m.  Thus 
again  we  see  that  the  index  of  damping  of  the  weight  function  is 
more  than  for  the  correlation  function,  and  only  at  large  X  do  these 
indexes  differ  little  from  one  another. 


t* 


A 

Fig.  5.  Change  of  e  depending  on  X  for 
wind  at  different  m. 


From  system  (28)  we  find  coefficients  C  and  D  for  optimum  inter¬ 
polation,  and  from  system  (29)  and  for  smoothing.  Simple 
formulas  are  obtained: 


C 


■  y  «  "Jv  »  U 

e  \€  —  e  1 


iV-O  ; 


C, 


(I r *  -  -  «*) 


(31) 

(32) 
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1  J 

Pig.  6.  Change  of  e  depending  on  A  for 
dew  point  at  different  m. 

III.  Solution  of  Integral  Equation  for  Weight  Function 

Let  us  turn  again  to  th/  case  when  the  amount  of  meteorological 
element  is  given  for  continuous  values  of  argument  x.  Finding  the 
weight  function  of  smoothing  f(y)  led  to  solving  integral  equation 
(14).  Among  various  means  of  solving  this  equation,  let  us  again 
examine  the  means  of  indeterminate  coefficients,  which  was  used  in 
the  second  section. 

We  set  again: 

r/r)  ~  e-«M,  r,(t)  — 

and  the  solution  we  look  for  in  the  form 

/W-^rw+AW, 

where  6(t)  -  designates  the  Dirac  delta- function. 


75 


Let  us  substitute  the  above  quantities  into  integral  equation 
(14).  After  this  it  is  possible  to  write: 


j  I*"' r’ I  IV"*  +  >M(y)l  dy  + 


+  j +  \A*e-*’  +  A,8(y)}  - 

M 


(33) 


Carrying  out  the  calculations,  we  compare  coefficients  of  e  , 
e  ,  e  .  We  obtain  a  system  for  determination  of  5,  A^,  A^. 


t 


0, 


(34) 

(35) 


From  system  of  equations  (34)  and  (35)  let  us  find  St  A.,  A. 

Z  t 

and  write  the  formulas  for  the  weight  function: 


c  i/TRrrnr 

.  . +nr-* 

i  _  o  _  2(«  —  j)  V 

1  (S  ->)={S  — S)- 


(36) 


We  set  (?E  to  be  small 


Let  us  examine  certain  limiting  cases. 

in  comparison  with  o^,  i.e.,  X  turns  into  zero.  Then  the  index  of 

damping  of  weight  function  S  coincides  with  the  index  of  camping  of 

2  2 

the  correlation  function  of  errors  0.  If  o  /a  grows  without  limit 

^  o 

(X  +  »),  then  S  approaches  a.  Obviously,  with  finite  X  we  have  the 
following  inequality:  a  <  S  <  0. 


Let  us  examine  now  the  particular  case  of  integral  equation 
(14),  when  errors  of  observation  bear  the  character  of  white  noise, 
and  the  upper  limit  has  a  finite  value,  rather  large  and  equal  to  N. 
Then  integral  equation  (14)  will  pass  into  a  Fredholm  equation  of 
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the  2nd  order.1 


fax  -  y)/(y)dy  +  >/(*)  -  rjx) . 


(37) 


We  assume  as  before  rg(T)  "  e-a^;  then  equation  (37)  is 
rewritten  in  the  form: 

f(y)dy+i.fix)-e-' 

6 

or 


]e-«*-»f{y)dy  +  (e f(y)dy +i./(x)-e-' .  ( 38 ) 

«  X 

The  solution  of  the  equations  (31)  we  reduce  to  the  solution  of 
a  differential  equation  with  constfant  coefficients.  For  this  let  us 
differentiate  twice  both -sides  of  (31)  with  respect  to  variable  x. 

We  write  the  result  of  the  first  and  second  derivatives: 


-« f  e-«*-»f(y)dy  +  a  [e«y*'f{y)dy + >.f(x\  «  -«  e-* . 

ii  (I 

**{  r-*-  *  f(y)dy  -  2*f(x) + >./'(*)  - 


(39) 


If  we  substitute  in  the  last  formula  (39)  the  magnitude  of  the 
integral  for  e-ouc  -  \f(x)  in  accordance  with  (38),  then  we  obtain 

).f(x)-,(2  +  l*)f(x)-0. 

whence  we  find  the  solution  of  the  differential  equation2 


‘Equation  (37)  can  be  reached  if  in  the  original  formula  (10) 
we  smooth  over  a  finite  segment  of  changes  x,  which  is  more  natural. 

In  other  words,  let  us  take  x)/{x)tnx).  On  the  otner  hand,  the 

assumptions  that  the  upper  limit  corresponds  to  the  fact  that  weight 
function  f(y)  turns  into  zero  when  y  i.  n. 

2Such  a  form  of  the  solution  was  accepted  by  us  for  the  discrete 
case  [see  formula  (18)]. 
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(40) 


where  S,=»]  Arbitrary  constants  L  and  M  will  be  defined 

below  from  boundary  conditions.1 


Let  us  designate 


j*~*r/(y)dy ■*  0  and  He,*/(y)dy—‘Q. 


set  in  formula  (38)  x  *  0  and  x  *  iif;  we  have  correspondingly: 


o>—  1. 

*  -  a’  fr/w,  -f  x/(iV)  —  r 


O  +  X/(0)-!t 

Qe-V +>./(, V)-e-*. 


Accurately  also,  setting  x  -  0  and  *  «  N  in  formula  (39),  for 
the  first  derivative  we  obtain  system  of  equations: 


-a  e-xQ  +  ).f(N)  -  -a  #-*. 


We  exclude  from  system  (41)  and  (42)  the  introduced  quantities 
Gt  Q,  we  write  the  system  for  the  boundary  conditions: 

,/(0, ->.,/«»-=  -2a.  , 


We  substitute  in  equations  (43)  the  values  of  f(0),  f( N)3  f'(0), 
and  f’(N)  from  formula  (40).  Then  we  go  to  a  system  of  equation 
for  determination  L  and  M: 

‘in  [8]  is  examined  a  more  general  case  of  the  solution  of  the 
Fredholm  equation  of  the  2nd  kind  with  a  symmetric  nucleus. 
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mtvmmm 


-<V  +  «)i-i-<7--  *)M  ^  —  x* 

*-•'<*-  7)I  +  -tf«.'(7  +  *)-0. 
>•  *  0. 


(4i*) 


Since  I  and  M  have  been  found  from  system  ( 44 ) ,  it  is  simple 
to  write  the  formula  for  the  weight  function 


/<y> 


__  2«(7  -r  .  »WT  — ») _ 

+  +  'k^<7  +  *)=  —  («— tH  * 


(45) 


Prom  this  formula  it  is  evident  that  f(S)  at  large  enough  N  approaches 
zero.  For  small  values  of  y  (at  large  N)  the  second  term  of  (45) 
is  smaller  than  the  first. 
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ON  THE  FILTRATION  OF  FIELDS  OBTAINED  AS  A  RESULT 
OF  NONLINEAR  TRANSFORMATIONS  OF  THE  INITIAL 
VALUES  OF  METEOROLOGICAL  FACTORS 


L.  N.  Strizhevskiy 


On  the  basis  of  the  theory  of  optimum  linear 
filtration  of  two-*dimensional  fields  we  calculate 
filtration  operators  for  fields  of  the  nonlinear 
part  of  vortex  advection  and  the  nonlinear  part 
of  the  individual  wind  derivative.  It  is  shown 
that  applying  results  of  calculations  to  a  simple 
scheme  of  wind  forecasting  obtains  stability 
before  an  unstable  system. 


Introduction 


As  is  known  original  meteorological  information  always  contains 
errors  connected  with  the  measurement  and  treatment  of  weather  data 
and  conditioned  by  the  superposition  of  a  whole  series  of  errors: 
errors  of  the  measuring  devices;  errors  of  the  observers;  errors 
which  appear  during  encoding,  decoding,  transmission  and  analysis  of 
meteorological  information  etc. 


These  errors  (for  analyzed  maps)  are  characterized  first  of  all 
by  a  radius  of  correlation  which  is  considerably  less  than  the 
correlation  radius  of  the  "true"  fields  of  meteorological  elements, 
which  allows  their  partial  division.  The  theoretical  basis  of  such 
a  division  is  the  theory  of  optimum  linear  filtration,  an  account 
of  which  can  be  found  in  the  survey  article  of  A.  M.  Yaglom  [7j. 
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For  the  first  time  this  theory  was  applied  to  a  meteorological 
problem  in  a  work  of  F.  Thompson  [4]..  In  [2]  the  operator  of 
optimum  smoothing  of  the  geopotential  $500  for  the  Pract*cal  corre¬ 
lation  function  of  the  observed  field  was  calculated.  This  work 
Indicated  that  during  linear  transformations  of  the  original  field 
of  the  optimum  filtration  operator  does  not  change. 

Otherwise  we  deal  with  the  acquisition  of  nonlinear  operators 
for  initial  fields.  Meanwhile  the  most  important  transformations  of 
meteorological  fields  used  in  the  process  of  numerical  forecasting 
and  weather  diagnosis  are  nonlinear. 

This  work  is  dedicated  to  finding  formulas  for  optimum  smoothing 
of  fields  obtained  as  a  result  of  certain  widespread  nonlinear 
transformations  of  original  fields. 

In  future  calculations  it  is  supposed  that  original  fields  are 
isotropic,  uniform,  distributed  normally,  and  the  mathematical 
expectations  are  equal  to  zero.  We  will  also  use  a  relationship 
given  in  [6],  enabling  us  to  express  the  moment  of  the  fourth  order 
of  a  normally  distributed  random  quantity  through  its  second  moments. 

J 

Filtration  of  a  Vortex  Advectlon  Field 


As  wa3  shown  in  [6],  the  spectral  density  of  the  nonlinear  part 
of  the  vortex  advectlon  can  be  recorded  as 


V*.  •  -  -sr/i  0  >*.  -  >•*.***  -  - 

- — 1|,  A*  —  v  ( 1 ) 


For  further  calculations  we  must  give  the  form  of  spectral 
density  functions  of  the  observed  geopotential  field  5^(1)  and 
noise  field  S  (A).  In  accordance  with  results  of  work  [6]  we  take 

JX 

-«•«» — y  , ,  (2) 

w  +  wf 


8  2 


2  -1 

where  the  parameters  take  values:  B^( 0)  »  230  dkm  ;  8  *  10  km 
Spectral  density  noise  we  give  as  In  [2]  in  the  form 


Setting  signal  and  noise  as  uncorrelated,  we  write 


Substituting  relationships  (2),  (4)  in  formula  (3)  and  intro¬ 
ducing  (used  in  [6])  a  change  of  variables: 


where 


i  “  nf-  + » «»(•  +  •)•  •y-4- »  *!«»{• 


and  the  polar  coordinates  in  plane  k: 

*i — *co««. 

We  obtain  the  expression  for  spectral  density  of  the  "true"  compo¬ 
nent  of  field  J(x,  y): 


o» 


*lB**CW*l> 


•  i'  •»* 

'.mtM  rr  — ->  Sf- 


■>»  <l«l 


jiB^CO*1*  ■ 
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The  first  integral  is  equal  to  the  right  side  of  formula  (1) 
and  expresses  the  spectral  density  of  the  field  J(x,  y)t  calculated 
according  to  observed  values  of  $(x,  y),  and  is  found  in  [6]: 

- 

l  «  (4p  + 

Finding  the  second  integral  also  is  not  difficult: 


d  •  *• 

j  sia*b  co*1  bd  b  j  d » -  8fl*  - *«•. 


The  third  and  fourth  integrals  in  expression  (6)  were  found 
numerically.  The  formula  for  calculating  the  smoothing  operator 
we  write  in  the  form 


where  JQ(Kp)  —  Bessel  function  of  the  1st  kind  of  zero  order. 

During  numerical  calculation  of  operator  K(p)  the  selection  of 
the  upper  boundary  of  integration  was  dictated  by  the  following 
considerations . 

v 

In  solving  forecast  problems  tfve  usual  assignment  of  meteoro¬ 
logical  fields  which  take  part  lr.  the  process  of  precomputation  in 
the  form  of  a  discrete  sequence  of  their  values  limits  from  above 
the  frequency  range  in  spectral  decompositions  of  these  fields.  In 
accordance  with  the  theorem  of  Kotel'nikov  [5]  maximum  wave  numbers 
for  a  right-angled  grid  in  the  x,  y  directions  are  determined  from 

*  « 

*w— 5-  *r~  S* 

where  h  t  h ^  are  the  grid  intervals  in  the  x,  y  directions.  For  a 
square  grid  it  is  apparent  that 


For  h  -  300  km  \ 


Results  of  calculations  according  to  formula  (7)  at  various 
values  of  parameters  a  and  B  (0).  characterizing  the  average  scale 
of  the  noise  field  and  the  mean-square  value  of  this  field  are  in 
Fig.  1. 


Fig.  1.  Weight  functions  for  smoothing 

the  nonlinear  part  of  vortex  advection 

B  (0)/BJO)  -  0.02  (a),  when  B  (0)/Bj0)  = 
n  ^  ti  tp 

*  0.1  (b).  1  —  a  •  3*10~3  km~x ,  2  —  a  * 

-  1.1-10'3  km'1,  3  -  a  -  0.7*10"3  km'1. 

In  order  to  use  the  obtained  function  K(p)  we  found  integrals 
of  form: 

•  * 


Replacing  values  of  the  neutralized  function  f(x,  y)  in  a  circle 
of  radius  and  in  rings  a^-a^  etc.  by  the  universal  means  in 

these  ranges  fjr,  y),— /,  y). ... .  we  obtain  the  working  formula 

for  smoothing  of  the  form 


y) 4- *«,-*/*.-*(-*.  y)+- 


(8) 


The  system  accepted  in  this  work  for  smoothing  is  seen  in  Pig.  2. 
The  following  working  formulas  for  smoothing  were  used: 


B„(0) 
B*(  0)  ' 
‘3,0  •  10- 


‘0,02 


KM' 


7  ■=  0,95/o-soo  +  0,OS7w-4»; 
a*=  I. MO-3  km-1, 

:..J  »*=  Q,877o~«).-r  0,  llJxn---™  -f-  OvO^/co-kjC'.v 
a  —  0,7  •  I0-3  km-1, 

7  <■=  0,657)- wn  -)■  0,2772«)-  «o  -f-  O.OSTuo-rat  • 


7 


II. 


a  «= 


fr,(0) 

fl*(0) 

3,0  •  10-* 


0,1 


KM“‘, 


/*=  0,9  l70-i»  +  0,09y»o-«9: 
a—>  i.l  •  10~a  km'1, 

=  0,797o-2oo  +  0, 197)00-450  +  0,027450-ho; 

a«0,7*  10~*  km"1, 

0,627o-mo  +  O,297)oo-450  -f  0,0974so-«o- 


J 

Pig.  2.  Partition  of  area 
during  the  smoothing  and 
numeration  of  points. 


Filtration  of  the  Nonlinear  Part  of  the  Individual 
Wind  Derivative 


Let  us  find  now  the  operator  of  optimum  smoothing  for  an 
expression  widespread  in  meteorological  applications 
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where  u,  v  -  components  of  wind  velocity . 


Let  us  present  F,  u3  v  in  the  form  of  stochastic  integrals: 

F(x,  y)-  kj, 

u(x,  y) —  f  dZfit,  X,), 

v(x,  y)  -  fj  g«^’>  dZjh,  ^).  ( 9 ) 


Differentiating,  we  write 

dZfa,  W\ ytZJtH.  **))• 

— -  m 


Substituting  the  components  of  wave  vector  v  with  the  aid  of 
relationships 


*,  —  —  A,  and  —  Af  —  kf. 

>  f 

Setting  it  and  v  as  independent  and  in  accordance  with  the  results 
of  [1]  spectral  densities  u  and  v  Identical,  we  obtain  the  following 
expression  for  spectral  density  F(x>  y): 


Sf  <*„  *,)  -  fj  l(».,  -  *,J*  +  (X,  -  *,)£(*,-  K  k2  -  WW.  ( 10 ) 


As  the  wind  correlation  function  we  took  an  expression  obtained 
in  [1] 


where 


BjiO)  —  250  j**b-*,  *  “  0,797, 

7-0,684,  /  —  1,28. 


(11) 


lFor  the  operator  F’(x,  y)  =  +  v^-,  obviously,  it  is  possible 


to  write 


k3)  m  j  j  K^i  —  k{jr  +  (2X,  —  *«)s)  S,  (X,,  Xf)  SM  (kt— kt,  kt^~X{fdXxdKt. 
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The  function  of  spectral  density  of  noises  was  assigned  in  the 

form 


$.<*> - (12) 

where  the  parameter  b  characterizes  the  position  of  maximum  spectral 
density. 

Considering  relationships  (4)  and  (10),  we  write  the  expression 
for  spectral  density  of  the  "true”  component  of  the  field  F(xt  y) 

Sr,(k„  *,)  -  J ]  [(2X,  -  ktf  +  *2  ~  WISJK  ^  (*,-X„ 

-« 

M 

-JJ  l(2X, -A,)»  +  (k1  -VJ  SJh.  kt^^d\d^- 

*—  «y 

-JJ  l(2x, -*,J*+(X.  -  WjiK  xjsfr.-x.,  *+ 

-w 

+  JJl(»i-*.J,  +  (X»-*,)*]Sld(Xlt  X^A.-X,,  A.-rX^X.dX*,  (13) 

where 

X,)  and  SJ}. ,,  A.,) 

are  determined  by  formulas  (11)  and  (12). 

The  first  integral  in  the  right  side  of  (13)  expresses  the 
spectral  density  of  field  F(x ,  y),  calculated  in  terms  of  observed 
values  of  u,  v. 

Obtained  as  a  result  of  numrrical  integration,  functions 

k„)  and  S„  (fe,,  k„)  were,  as  one  would  expect,  even.  Con- 
F  1  2  t  g  1  2 

sidering  this,  and  likewise  the  anisotropism  of  these  functions,  let 
us  write  down  the  expression  for  optimum  smoothing  of  field  F 

K(x,  y)«4irjj  cos(M 4- kjy)  dkxdkt.  (14) 
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The  contour  lines  of  function  K(xt  y)  for  some  values  of  a  and 
b  are  given  in  Pig.  3  (the  picture  obviously  Is  symmetric  with 
respect  to  the  axes  of  coordinates). 


Pig.  3.  Weight  functions  K(x ,  y)  for  the 

smoothing  of  field  F  when  B  (0)/B  (0)  • 

n  u 

-  0.04  and  a  -  0.7*10"3  km"1. 

> 

The  most  interesting  result  in  this  case  is  an  anisotropism  of 
the  smoothing  operator,  which  is  explained  by  the  different  character 
of  change  in  spectral  densities  of  signal  and  noise  (different  degree 
of  anisotropism)  during  the  transformation  (12).  We  must  note  that 
from  the  anisotropism  of  conversion  (10)  follows  immutable  aniso¬ 
tropism  of  K(x,  y):  any  linear  anisotropic  conversion,  for  example, 
differentiation  according  to  one  of  the  variables,  retains  the 
isotropism  of  the  smoothing  operator. 

The  degree  of  averaging  in  the  direction  of  the  A'  axis  with 
certain  characteristics  of  noises  is  almost  one  and  a  half  times 
greater  than  the  degree  of  averaging  along  the  7  axis  near  the  point 
to  which  the  value  of  the  neutralized  field  belongs.  With  distance 
from  the  origin  of  coordinates  function  K(xt  y)  becomes  more 
isotropic. 

Smoothing  formulas  were  obtained  by  means  analogous  to  those 
described  in  the  previous  section.  Let  us  give  the  results  of 
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calculations  for  certain  values  of  noise  parameters  (the  arrangement 
of  points  is  given  in  Pig.  1): 


fl.(O)  '• 

I.  4-2.5-I0-*  KM-' 

1)  0=0,7. 10-3  KM~\ 

7-  0,72/\>  +  0.06^,  +  F>)  +  0.02(f,  +  f,  +  F>  +  F,  + 

+  F,  +  F.+  Fr  +  M: 

J)  u-0.4  10-*  km-'. 

F-  0.77 Ft  •+■  QfloiFy  t  Ft)  +  0.016 {Ft  +  F>  +  F„  +  Ft+  . 

F-.  +  F\  +  Fr  +  Fr). 

II.  b  -2,25- 10-*  km-' 

1)  a  —  0.7-10-3  km-1, 

7-.0,5lf0  +  0.  !25</r,  +  /•’,)  +  0.022 (f,  +  F,  +  F\  +  Ft +F;  +  F,+ 

+  Fr  +  F «•)  -j-  0.0l(P,.  +  Fr  4-  Fy  -(-  Fr  +  Fr  4*  Fr)\ 

2)  a  —  0.4-10-*  km-1, 

F- 0,55  F*  4-  0, 1  \(Ft  4-  Ft)  4-  0,02 (f ,  +  Ft  + Fk  + Ft  + F,  +  Ft  + 

+  Fy  4-  Fr)  4-  0.0l<f ,  •  4-  4-  Fy  4-  Fr + /V  4-  Frb  (15) 


Smoothing  formulas  (15)  show  that  the  degree  of  smoothing 
depends  mainly  on  the  magnitude  of  parameter  b  or,  in  other  words , 
on  the  degree  of  overlap  of  the  spectras  of  signal  and  noise.  With 
an  increase  in  the  overlap  naturally,  it  increases. 


Using  Filtration  in  the  System  of  Wind  Forecastinc 


Results  obtained  in  the  previous  section  were  used  in  the  one- 
level  system  of  wind  forecasting  described  in  [3]. 


Here  briefly  is  its  variant,  using  as  raw  data  the  fields  of 
actual  wind. 


Motion  equations  are  written  in  the  form: 


4* 

41 


(16) 


where  u',  v'  are  the  components  of  the  deviation  of  wind  from 
geostrophic ; 
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to  _  to  1  „  to  ,to 
to  _  to  ,  to  ,  to 

Tr“-jr+*Tjr+  TT 


The  remaining  designations  are  generally  accepted.  The  Ageo- 
strophic  part  of  the  stream  function  is  connected  to  v'  by  the 
relationships : 

•'-4^  UT) 


Values  of  can  be  determined  according  to  the  field  of  prac¬ 
tical  wind  from  the  balance  equatior 

.  a./  _  /  to  \* ,  /  to  \*  .  «  to  to  _  to  ,  _  to 


,  , ,  /to  ,  /  to  \*  ,  0  to  to  ■  to  ,  _  to 


Further,  with  the  aid  of  relationships  (17)  we  find  u'  and  v' 
and  integrate  equations  (16)  o^tr  time.  The  integration  was  in 
Lagrangian  coordinates,  which  was  connected  with  interpolation  of 
the  values  of  wind  components  every  grid  interval  in  the  inside 
areas  of  the  grid  squares.  The  latter  led  to  noticeable  leveling 
of  results. 

Attempts  to  integrate  these  equations  in  Euler  coordinates  with 
the  aid  of  relationships: 


-  «<«’  +  « t(lv,w-  -  *•»  -^J. 

*'»  —  tr<*>  -  Itf  -f  V'** 


where  n  is  the  number  of  the  time  interval,  were  unsuccessful  -  in 
the  absence  of  smoothing  the  system  was  unstable.  Its  stability  was 
achieved  after  a  number  of  experiments  by  smoothing  in  accordance 
with  the  above  results.  On  every  time  step  the  smoothing  was  carried 
out  twice:  the  sums  of  the  nonlinear  terms  in  the  right  side  of 
equations  (18)  by  formula  (15),  while  the  most  successful  results 
were  obtained  with  the  following  values  of  parameters: 
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0,04,  *  -  2.5  •  I0-*  km-‘. 


and  weak  smoothing  of  wind  components  according  to  formula 

«- 0.94* +  0.015  ija,. 


where  u  is  the  wind  component . 


Results  of  the  estimation  of  forecast  quality  are  in  Table  1, 
where  |C*  -<£|  -  mean  absolute  error  of  wind  prognostication,  !aC„|  - 
mean  modulus  of  vectorial  wind  variability  for  twenty-four  hours. 


Integration  of  forecast  equations  for  formulas  ( 1 8 )  using 
precalculated  smoothing  somewhat  improved  results  in  comparison  with 
the  earlier  variant  systems,  which  is  connected  in  our  opinion 
basically  with  obtaining  more  practical  values  of  the  absolute 
values  of  wind  velocity.  Furthermore,  utilization  of  formulas  ( 1 8 ) 
considerably  simplified  the  computer  program  and  shortened  calcu¬ 
lation  time  in  comparison  with  the  previous  variant. 


Table  1. _ 

;  Bilinear  interpolation  Optiaur  nitration 


)  (Langrangian  coordinates)  | 

(Euler  coordlnatea) 

Data  | 

V 

1 

K-* 

e.-ej 

(  ■/•  ) 

! 

|A?*i 

(  a/e  ) 

Laval  700  abar 


0300  2  Mar  I960 

6.7 

OM 

M 

0M 

0300  3  fUr  1960 

M 

0,76 

U 

0.7J 

C300  4  Her  I960 

7.0 

a  77 

a 

m 

0300  5  Mar  1960 

74 

an 

u 

w* 

0300  26  Fab  1961 

U 

043 

y 

0M 

1500  29  Fab  1961 

IM 

0,75 

AXO 

47« 

Laval  750  abar 


0300  3  Mar  1960 

X* 

043 

u 

04* 

0300  4  Hu-  1960 

64 

1.10 

M 

144 

0300  5  ffer  I960 

3.7 

0.73 

X2 

tm 

0300  26  t%r  1961 

U 

048 

to 

Of! 

1500  26  Fab  1961 

u 

047 

M 

046 

Average . 

040 

04> 

U3 

477 
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OBJECTIVE  ANALYSIS  OP  THE  QEOPOTENTIAL  FIELD  USING 
SUPPLEMENTAL  INFORMATION 

V.  D.  Sovetova 

Two  methods  of  objective  analysis  are 
described.  The  first  method  is  dedicated  to  the 
analysis  of  the  isobaric  surfaces  of  the  tropo¬ 
sphere  over  the  northern  hemisphere,  and  the 
second  to  analysis  of  stratospheric  levels  over 
Eurasia.  To  improve  the  quality  of  analysis 
over  incompletely  treated  regions  and  at  high 
altitudes  additional  information  is  used. 

> 

Numerical  analysis  and  forecasting  the  baric  field  over  the 
northern  hemisphere  is  hampered  in  regions  incompletely  covered  by 
meteorological  information.  The  quality  then  drops  off  [7]  over 
water  surfaces  or  the  Pacific,  Atlantic,  Indian  Oceans,  and  over 
Africa,  where  there  are  few  points  of  radiosounding  of  the  atmo¬ 
sphere.  On  the  other  hand,  for  example,  the  territory  of  Eurasia 
on  the  whole  is  well  covered  by  the  data  of  observations  in  the 
troposphere  and  unsatisfactorily  in  the  stratosphere.  Therefore  for 
improvement  of  analysis  in  such  cases  it  is  expedient  to  use  addi¬ 
tional  information.  Two  methods  of  restoration  and  analysis  of  the 
geopotential  field  are  considered:  the  flrts  is  developed  for  the 
troposphere,  and  the  second  for  the  strotosphere . 

The  first  problem  is  solved  for  the  northern  hemisphere  (grid 
of  1545  points,  intervals  at  470  km)  and  for  four  levels  of  the 
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troposphere  -  1000,  850,  500  and  300  mbar.  At  the  surface  of  the 
earth  there  were  sufficient  data  for  a  successful  analysis  In  almost 
all  regions  of  the  northern  hemisphere.  The  quantity  of  information 
sharply  falls  upon  transition  to  the  higher-lying  isobaric  surfaces 
(850,  500,  300  mbar),  where  there  are  data  only  for  the  points  of 
atmospheric  sounding.  To  supplement  information  at  these  levels  one 
could  use  the  data  of  measurements  of  the  surface  field  of  pressure. 
Por  this  purpose  for  interpolation  onto  points  of  a  regular  grid, 
except  radiosonde  stations  (n  -  589),  station  of  making  only  ground 
observations  (n  •  290)  have  been  used.  Over  the  oceans  sources  of 
additional  information  are  observations  of  ^8  ships  of  the  merchant 
marine . 

The  quantity  of  information  on  100  and  850  mbar  was  substantially 
different,  therefore  before^objective  analysis  were  reconstructed 
absent  measurements  of  the  geopotential  at  850  mbar  according  to 
available  data  of  surface  pressure;  the  absent  information  for 
300  mbar  was  reconstructed  according  to  data  for  500  mbar. 

Reconstruction  of  unavailable  information  was  conducted  by  the 
method  of  spatial  optimum  interpolation  [3]  using  mutual  correlation 
functions  of  the  geopotential  [8]  for  1000  and  850  mbar,  500  and 
300  mbar.  Reconstruction  could  have  used  simultaneously  the  data  of 
the  lower  and  original  levels.  However,  because  information  at  the 
original  level  as  a  rule  is  absent  in  those  points  around  which  at 
a  distance  of  the  radius  of  correlation  there  are  no  data  of  obser¬ 
vations,  for  its  reconstruction  the  information  of  five  stations  of 
one  lower  level  better  covered  by  meteorological  information  was 
used.  One  of  these  five  stations  had  to  lie  on  one  vertical  line 
with  the  station  of  the  original  level,  and  the  remaining  four  were 
selected  such  that  each  of  them  covered  one  of  the  four  quadrants 
arranged  around  the  reference  station.  When  in  the  circle  wi*-h  the 
correlation  radius  it  was  not  possible  to  find  four  symmetric 
stations  in  which  there  would  be  data  of  measurements,  any  four 
stations  were  selected,  independently  of  their  arrangement  around  a 
point.  But  even  if  this  condition  could  not  be  made,  reconstruction 
used  the  information  of  only  one  lower  point. 
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As  a  result  of  the  reconstruction,  subsequently  In  the  inter¬ 
polation  onto  regular  grid  points  of  the  geopotential  of  1000  and 
850  mbar  879  stations  can  take  part;  objective  analysis  of  the  baric 
field  of  500  and  300  mbar  is  carried  out  from  the  data  of  589 
radiosonde  stations. 

To  estimate  the  quality  of  the  reconstruction  by  the  above 
means  a  numerical  experiment  was  conducted  by  means  of  the  artificial 
exclusion  of  the  entire  information  on  850  mbar.  The  mean  square 
error  of  reconstruction  of  the  geopotential  of  850  mbar  from  the 
data  of  0300  21  January  1964  was  4.6  dam  over  all  the  northern 
hemisphere,  3-3.  dam  over  North  America, and  4.3  dam  over  western 
Europe.  The  geopotential  of  500  mbar  was  reconstructed  analogously 
from  data  for  850  mbar  with  a  mean  square  error  of  5.0  dam;  data  for 
300  mbar  was  reconstructed  from  data  for  500  mbar  with  a  mean  square 
error  of  7.6  dam.  Taking  into  account  that  under  practical  condi¬ 
tions  reconstruction  of  information  is  necessary  only  in  the  inade¬ 
quately  covered  regions  and  the  obtained  errors  are  somewhat  less 
than  the  corresponding  deviations  from  climatic  standards  which  must 
be  used  here,  the  results  of  the  reconstruction  can  be  considered 
satisfactory. 

Comparison  of  standardized  autocorrelation  functions  of  the 
geopotential  showed  the  proximity  of  functions  for  1000  and  850  mbar 
and  also  the  proximity  of  the  functions  of  50C  and  300  mbar.  The 
difference  of  correlation  functions  between  the  first  and  second 
pairs  of  surfaces  is  rather  substantial.  In  connection  with  this 
the  objective  analysis  of  the  reconstructed  fields  was  carried  out 
in  two  stages:  is  from  the  beginning  for  lsobaric  surfaces  of  1000 
and  850  mbar,  after  that  for  another  pair  of  surfaces  -  500  and  300 
mbar. 


The  objective  analysis  was  conducted  by  the  method  of  optimum 
interpolation  [4],  Not  describing  all  units  of  the  program,  let  us 
stop  only  on  its  separate  parts,  determining  the  distinction  of  our 
methodologies  from  the  methodology  accepted  at  the  Hydrometeoro¬ 
logical  Center  of  the  USSR  [6]. 


1.  For  the  characteristic  features  of  the  average  climatic 
fields  of  the  geopotential  (of  standards)  on  1000  and  500  mbar  were 
assigned  fields,  calculated  in  NIIAK  [1].  Standards  for  two  other 
surfaces  (850  and  300  mbar)  were  calculated  by  the  formula: 

■Zm*  “  £|«m  +  A  Zmk 

where  “Kz  -  average  climatic  value  of  the  thickness  of  the  layer 
between  the  corresponding  isobarlc  surfaces,  figured  according  to 
the  formulas: 

for  winter 

SI mm  -  -  32.5<  *ln  *)*  +  7,5  st»  *  +  140. 

IzS~-48,0i*ln*y*  +  3«>; 

\ 

for  summer 

AZSL -  -$8.5<»ln -f  27.5 tin  f  +  140, 

SIS  -  -3*.5<«ln  rf  +  7,5  tin  ?  +  385. 

Comparison  of  actual  standards  of  the  geopotential  for  850  and 
300  mbar  with  calculated  standards  showed  their  close  correspondence. 

2.  The  stations  were  situated  in  such  a  way  that  their  ordinal 
numbers  increased  along  the  band  and  from  one  band  to  another.  The 
first  station  of  every  band  took  on  the  feature  by  which  the  bands 
were  estimated.  During  interpolation  onto  the  point  with  coordi¬ 
nates  x,  y  the  number  of  bands  was  calculated  beginning  from  the 
first,  and  the  sum  was  equated  with  y  -  3.  The  comparison  determined 
in  which  half  of  the  bands  (right  or  left)  the  point  is  located. 

Only  those  stations  which  were  located  in  the  same  half  as  the  point 
and  which  were  located  in  the  band  between  (y  -  i)  and  (y  *  3 )  were 
tested  to  see  if  they  belonged  to  the  neighborhood  of  the  considered 
point.  As  soon  as  they  began  to  follow  stations  with  coordinates 
more  than  (y  *  3),  the  search  ceased. 


3.  The  interpolation  onto  regular  grid  points  was  conducted 
in  the  following  manner,  ^rom  the  beginning  an  attempt  was  made  to 
select  the  eight  nearest  stations,  two  stations  in  each  quadrant.  If 
this  requirement  was  not  fulfilled,  eight  asymmetrical  stations  were 
used.  If  in  the  whole  square  of  6  *  6  grid  intervals  there  were  not 
eight  stations  with  data  from  observations,  the  choice  went  to  four 
stations  satisfying  the  requirement  of  symmetricity .  If  four 
asymmetrical  stations  were  lacking,  only  two  stations  were  used.  To 
calculate  the  weight  coefficients  a  system  of  equations  of  the  eighth, 
fourth  or  second  order  was  made  up,  depending  on  how  many  stations 
with  measurement  data  were  found. 

On  the  edges  of  the  map  frequently  Iz  was  not  possible  to  find 
in  a  given  square  even  two  stations.  In  this  instance  in  points  on 
1000  and  850  mbar  the  corresponding  climatic  values  were  used. 

Usually  these  points  were  in  the  southern  latitudes,  where  the 
dispersions  were  small,  and  that  is  why  the  absolute  errors  of  calcu¬ 
lation  were  small. 

On  500  and  300  mbar  the  interpolation  onto  points  was  done  from 
eight  or  four  stations.  With  fewer  stations  interpolation  was 
replaced  by  extrapolation  from  the  lower  surface  analogously  the 
reconstruction  of  lacking  information.  The  only  difference  is  that 
In  the  first  case  the  information  was  reconstructed  for  a  station, 
and  In  the  second  case  for  a  regular  grid  point  according  to  data 
for  a  point  located  on  the  underlying  surface  (850  mbar).  After 
this  procedure  the  value  of  the  geopotential  in  a  point  at  300  mbar 
was  calculated  by  means  of  extrapolation  to  the  upper  level  of  the 
value  in  the  same  point  of  500  mbar.  In  these  cases  spatial  extrapo¬ 
lation  of  deviations  from  standard  along  one  lower  point  was  con¬ 
ducted. 

Comparison  of  results  of  objective  and  synoptic  analyses 
obtained  mean  square  errors  of  comparison  for  two  types  of  points. 

The  first  type  includes  points  to  which  Interpolation  was  produced 
according  to  the  data  of  the  considered  level;  the  second  type 
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includes  points  to  the  values  of  geopotentials  which  showed  the 
effect  of  extrapolation  from  below  or  these  values  are  equal  to  the 
middle  climatic  value.  Errors  in  decameters  are  presented  in 
Table  1. 


Table  1. 

Surface,  mbar .  1000  850  500  300 

a  dam: 

first  type .  1.7  2.8  2.9  6.5 

second  type .  2.1  2.9  6.0  8.9 


As  can  be  seen  from  the  tables,  using  standards  on  the  edges 
of  the  map  instead  of  interpolated  values  increases  the  eiror  of 
comparison  by  0.4  dam  on  1000  mbar;  the  increase  of  error  (by  0.1 
dam)  on  850  mbar  is  entirely  insignificant.  On  500  and  300  mbar  the 
error  of  comparison  substantially  increases  for  cases  of  extrapol¬ 
ation  from  below  (to  3.1  and  2.4  dam,  respectively).  This  increase 

/ 

can  be  explained  uniquely,  because  it  is  conditioned  not  only  by  the 
difference  of  methods  of  calculation,  but  also  by  the  considerable 
drop  in  quality  of  synoptic  analysis  at  these  levels. 

Comparison  of  the  amounts  of  deviations  (Table  1)  with  standard 
deviation,  which  is  given  in  a  work  of  S.  A.  Mashkovich  and  S.  I. 
Gubanova  [7],  indicates  the  better  quality  of  objective  analysis  of 
1000  mbar  geopotential  surfaces  for  our  method.  Here,  obviously, 
we  see  the  effect  both  the  using  the  additional  information,  and 
using  the  autocorrelation  function  of  the  considered  level  in  calcu¬ 
lations  instead  of  the  function  of  the  geopotential  at  500  mbar  for 
all  levels. 

To  estimate  the  effect  on  numerical  analysis  of  the  additional 
information  over  sea  water  surfaces  the  following  experiment  was 
conducted.  The  data  of  ships  of  the  merchant  marine  were  artifi¬ 
cially  excluded  and  the  obtained  analysis  was  equated  with  the 
analysis  of  the  usual  variant.  Mean  square  errors  0,  obtained  during 
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the  comparison  of  these  two  analyses,  are  given  in  Table  2.  Here 
is  shown  the  number  of  points  n  in  which  was  discovered  the  differ¬ 
ence  in  these  analyses. 


Table  2. 


Surface,  mbar.... 

1000 

850 

500 

300 

a  dam .  . 

. . . .  9.2 

9.2 

10.6 

12.9 

n . . .  . 

....  214 

207 

211 

211 

As  can  be  seen  from  this  table,  measurements  of  pressures  on 
48  ships  show  the  effect  on  values  of  the  geopotential  in  large 
quantities  of  regular  grid  points,  and  the  difference  in  analyses  is 
rather  considerable.  For  1000  mbar,  where  the  oceans  are  covered 
considerably  better  by  the  information  than  at  heights  and,  conse¬ 
quently,  synoptic  analysis  is  more  reliable,  the  objective  analysis 
in  the  absence  of  data  from  additional  ships  was  compared  with 
synoptic  analysis.  It  turned  out  that  the  absence  of  additional 
ships  increases  the  mean  square  error  of  objective  analysis  over 
oceans  in  comparison  with  synoptic  analysis  by  8.0  dam,  i.e.,  it  can 
be  considered  that  the  effect  of  additional  information  on  the 
results  of  numerical  analysis  is  rather  considerable. 

During  the  analysis  of  calculated  fields  of  pressure  over  all 
the  northern  hemisphere  it  was  explained  that  in  the  Inadequately 
covered  water  surfaces  of  the  oceans  during  extrapolation  from  lower* 
levels  to  500  and  300  mbar  deeper  cyclones  than  noted  on  the  oper¬ 
ational  constant-pressure  charts  are  obtained.  The  authenticity  of 
the  calculated  values  Is  difficult  to  verify,  because  synoptic 
analysis  in  these  regions  is  accomplished  according  to  the  data  of 
the  observations  of  one  or  two  stations.  However,  the  wind  data 
allow  proposing  the  presence  here  of  gradients  greater  than  noted 
on  the  map.  For  example,  21  January  1964  in  the  center  of  a  cyclone 
over  the  Atlantic  Ocean  (42°  north  lat.  and  37°  west  long.)  our 
calculations  obtain  a  value  of  the  geopotential  for  500  mbar  of 
512  dam.  On  a  constant-pressure  map  prepared  by  the  forecasters  at 
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OMTs  SS3R  [Hydrometeorological  Center  USSR]  in  the  center  of  this 
cyclone  a  height  of  517  dam  ia  noted;  on  the  map  of  another  fore¬ 
casting  institution  in  the  center  of  the  cyclone  a  height  of  530  dam 
is  noted.  Thus,  our  data  differ  in  the  first  case  by  5  dam,  and  in 
the  second  by  18  dam.  The  actual  depth  of  the  cyclone  and  its 
accurate  location  are  unknown  inasmuch  as  in  the  range  of  the  cyclone 
there  are  observations  of  only  individual  stations. 

The  second  problem  was  solved  for  levels  of  the  stratosphere 
over  the  territory  of  Europe  and  Asia  (29  *  22  grid,  interval  ^50  km) 
The  primary  surface  whose  data  (geopotential  Z  and  temperature  t) 
were  used  for  the  reconstruction  of  information  on  higher  levels  was 
200  mbar.  Selecting  this  surface  avoids  the  effect  of  tropopause; 
losses  of  information  at  this  level  in  comparison  with  300  mbar  are 
insignificant  (2-3$). 

Analysis  of  our  autocorrelation  functions  of  the  geopotential 
for  stratospheric  levels  showed  a  substantial  difference  between 
functions  for  100  mbar  and  for  300  and  200  mbar.  Therefore  analysis 
at  100  mbar  was  conducted  using  the  autocorrelation  functions  corre¬ 
sponding  to  this  level. 

After  monitoring  the  original  Z  and  t  data  [5],  the  missing  Z 
and  t  information  was  reconstructed  at  the  primary  200  mbar  level  by 
the  method  of  optimum  interpolation  [4]  according  to  the  observations 
of  eight  stations.  Values  at  stations  located  on  the  edges  of  the 
maps  were  reconstructed  in  some  other  way.  For  these  stations  (n  = 

=  60)  separately  from  the  general  field  the  field  of  pressure  and 
temperature  was  assigned  according  to  a  previous  map  (radiosonde 
data  if  available,  or  values  taken  from  isohypses).  In  the  absence 
of  information  after  a  similar  period,  the  values  of  pressure  and 
temperatures  at  the  regional  stations  their  values  in  a  previous 
period  of  observations  were  taken.  As  a  result  of  reconstruction 
on  200  mbar  we  had  data  about  geopotential  and  temperature  in  ^39 
points . 
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Using  further  the  mutual  correlation  functions  of  temperature 
[2],  by  the  method  of  spatial  interpolation  [3]  the  missing  tempera¬ 
ture  was  reconstructed  on  100  mbar  (note  that  the  missing  information 
at  this  level  comprises  15-20?  of  the  relative  surface  at  300  mbar, 
not  including  erroneous  data).  Then  according  to  the  average 
temperature  of  the  layer  the  height  of  the  surface  at  100  mbar  was 
resolved  in  the  points  where  there  was  no  or  rejected  information. 

Before  objective  analysis  all  available  information  was  again 
monitored.  Rejected  data  were  reconstructed  by  the  method  of  optimum 
interpolation  on  a  fixed  level.  Because  the  accepted  methodology 
of  reconstruction  of  information  allows  having  at  all  levels  always 
the  same  number  of  data,  it  is  possible  to  attach  at  every  regular 
grid  point  the  very  best  located  stations  according  to  the  data  of 
which  the  interpolation  will  be  carried  out.  In  the  first  place 
objective  analysis  was  conducted  for  the  central  rectangle  (24  x  16 
points),  then  for  surrounding  rectangles  using  already  calculated 
values  in  certain  points  lying  on  the  contour  of  the  central  rec¬ 
tangle. 
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Improving  the  Methodology  of  Forecasting  Baric 
Field  to  Several  Days 

Mashkovich  S.  A.,  Transactions  GMTs,  1968, 

No.  19,  pp.  3-9 

A  generalization  is  made  of  a  prediction  system  proposed  earlier 
by  S.  A.  Mashkovich  and  Ya.  M.  Kheyfets  (Transactions  TsIP,  No.  93). 
The  generalization  consists  of  recording  the  effect  of  horizontal 
turbulent  mixing  in  the  vortex  equation  and  the  effect  of  surface 
friction.  The  latter  is  introduced  when  the  boundary  condition  at 
sea  level  for  the  vertical  component  of  velocity  is  recorded.  The 
original  equations  are  vortex  and  beat  influx  equations,  linearized 
relative  to  the  zone  of  motion.  A  three-level  model  of  the  atmo¬ 
sphere  is  considered.  Motion  is  assumed  to  be  quasi-solenoidal.  The 
solution  for  the  stream  function  is  looked  for  in  the  form  of  a  series 
in  spherical  functions.  The  problem  was  solved  by  computer.  Seven¬ 
teen  series  of  forecast  maps  for  a  period  up  to  four  days  were  cal¬ 
culated.  The  quality  of  the  forecasts  for  four  days  is  estimated. 
Thanks  to  the  generalization  of  the  system  it  was  possible  to 
diminish  the  absolute  forecast  error  on  the  average  by  15-201.  In 
July  1966,  the  composition  of  forecasts  of  ground  pressure  by  this 
system  became  operational. 

Table  1,  Ill.  1,  Bibl.  3. 


One  Finite-Difference  Algorithm  for  the  Solution  of 
the  Vortex  Equation  for  the  Middle  Troposphere 
Over  the  Northern  Hemisphere 

Isayev  N.  V.  and  Fuks-Rabinovich  M.  S.,  Transactions 
GMTs,  1968,  No.  19,  pp-  10-21 

A  system  for  forecasting  the  geopotential  at  the  middle  tropo¬ 
sphere  over  the  northern  hemisphere  is  examined.  To  prevent  non¬ 
linear  instability  a  finite-difference  approximation  of  nonlinear 
terms  was  used,  proposed  by  A.  Arakava.  With  such  an  approximation 
there  is  conservation  of  the  quadratic  integral  features,  i.e., 
conservation  on  the  forecast  range  of  the  integrals  of  the  vorticity 
rate,  its  square  and  kinetic  energy.  Integration  in  time  was  con¬ 
ducted  using  several  methods,  among  which  most  suitable  was  the 
Adams  method. 

The  system  was  used  to  execute  calculations  over  a  long  period 
In  order  to  steady  the  character  of  change  in  kinetic  energy  of 
forecast  fields  In  time.  It  turned  out  that  the  finite- 

difference  algorithm  allows  calculation  of  a  forecast  over  long 
periods,  moreoveT  the  kinetic  energy  of  forecast  fields  during  cal¬ 
culation  remains  practically  constant.  Results  of  testing  the 
system  over  a  period  of  up  to  three  days  are  presented,  and  a 
quantitative  and  qualitative  analysis  of  the  obtained  forecast  fields 
are  given. 

Table  3,  Ill-  5,  Bibl.  14. 
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Experiment  In  the  Numerical  Solution  of  Balance  Equations 
Within  the  framework  a  Quaal-Solenoidal  System  of 
forecasting  the  Oeopotentlal  on  the 
Northern  Hemisphere 

Sitnikov  I.  0.  and  Krichak  S.  0.,  Transactions 
GMTs,  1968,  No.  19,  pp.  22-30 

A  finite-difference  system  for  solving  the  balance  equation  on 
the  northern  hemisphere  Is  given,  and  a  number  of  characteristic 
features  found  with  the  solution  of  this  equation  within  the  frame¬ 
work  of  the  quasi-solenoidal  system  of  forecasting  the  geopotential 
at  the  middle  troposphere. 

For  solving  the  balance  equation  the  method  of  successive 
approximations  is  used,  and  the  balance  equation  is  reduced  to  a 
Poisson  equation  relative  to  the  stream  function  and  values  of  the 
right  side  are  determined  from  the  previous  approximation.  One 
feature  of  the  solution  of  the  balance  equation  is  the  reduction  of 
the  rate  of  convergence  of  the  Iteration  process  in  the  low  latitudes. 

Obtained  as  a  result  of  solving  the  balance  equation,  the  field 
of  the  stream  function  is  used  as  the  original  for  the  quasi- 
solenoidal  forecast  based  on  a  baratropic  model.  At  the  end  of  each 
forecast  period,  i.e.,  after  24,  48  and  72  hours,  the  geopotential 
field  is  located  by  inversion  of  the  balance  equation,  i.e.,  solving 
it  relative  to  the  geopctential  according  to  the  known  distribution 
of  the  stream  function.  f 

An  example  of  calculating  a  forecast  from  the  quasi-3olenoidal 
model  is  presented,  and  it  is  compared  with  a  quasi-geostrophic 
forecast. 

Ill.  6,  Blbl .  9. 


About  Research  on  Warming  Trends  In  the  Stratosphere 
Using  Numerical  Experiments 

Barg  B,,  Hashkovich  S.  A.,  Transactions  GMTs, 

1968,  No.  19,  pp.  31-43 

The  evolution  of  a  narrow  band  of  temperature  perturbation  in 
the  stratosphere  Is  studied.  The  starting  equations  are  those 
proposed  by  A.  S.  Dubov  (Transactions  GGO,  No.  124).  The  Investi¬ 
gation  is  conducted  according  to  stylized  Initial  conditions:  it 
assumes  that  the  original  temperature  perturbation  has  a  circular 
form  and  occupies  a  limited  range.  An  initial  disturbance  of  the 
baric  field  is  absent.  Finite-difference  approximation  of  equations 
is  described,  a  numerical  method  of  solution  is  formulated,  the 
calculating  stability  of  the  solution  is  investigated. 
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Preliminary  results  are  given  on  the  basis  of  which  the  following 
basic  conclusions  are  made: 

1)  the  range  of  temperature  perturbation  used  in  the  calcula¬ 
tions  does  not  substantially  affect  circulation  in  the  stratosphere, 
specifically,  the  polar  vortex  does  not  break  down; 

2)  perturbation  in  the  temperature  field  also  causes  a  pertur¬ 
bation  of  the  baric  field,  the  range  of  which  grows  in  the  consid¬ 
ered  period  of  time  (4-6  days),  while  the  range  of  the  temperature 
perturbation  diminishes; 

3)  maximum  pressure  perturbation  remains  during  the  first  48 
hours  at  the  same  latitude  on  which  the  temperature  perturbation  lay, 
then  the  geopotential  perturbation  shifts  to  the  north,  while  the 
temperature  perturbation  remains  at  the  original  latitude. 

Table  1,  Ill.  4,  Bibl.  10. 


Application  of  the  Theory  of  Filtration  of  Random 
Processes  to  Certain  Problems  in 
Objective  Analysis 

Ovchinskly  B.  V. ,  Transactions  GMTs,  1968, 

No.  19,  pp.  44-57 

Formulas  are  set  for  calculating  weight  functions  of  optimum 
interpolation  and  smoothing  meteorological  elements  when  correlation 
function  of  the  true  magnitude  (signal)  and  observation  errors 
(noise)  are  known. 

During  optimum  Interpolation  one  must  take  Into  account  that 
the  data  of  observations  contain  random  errors,  and  an  attempt  Is 
made  to  achieve  a  magnitude  of  the  meteorological  element  in  a 
regular  grid  point  free  from  these  random  errors.  If  the  examination 
Includes  small-scale  fluctuations,  then  it  Is  necessary  to  reject 
the  assumption  of  noncorrelatabllity  of  the  random  errors.  This 
circumstance  was  pointed  out  by  Thompson,  who  used  the  theory  of 
filtration  of  random  processes  for  optimum  smoothing  of  meteoro¬ 
logical  fields. 

The  structure  and  features  of  weight  functions  were  studied 
both  for  discrete  and  continuous  data  of  observations.  The  problem 
of  determining  weight  functions  reduces  to  solving  the  Integral 
Wlener-Hopf  equation  (or  a  discrete  model  of  the  Wiener-Hopf  equation 
tlon).  Of  the  various  ways  of  jolvlng  this  equation  the  method  of 
Indeterminate  coefficients  was  used;  in  this  case  the  correlation 
function  is  approximated  by  the  sum  of  the  model  functions. 

Cases  of  a  solution  when  errors  jf  observation  are  "white  noise" 
or  when  the  number  of  stations  increases  without  limit  are  examined 
also . 


Ill.  6,  Bibl.  11. 
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Filtration  of  Fields  Obtained  as  a  Result  of  Nonlinear 
Transformations  of  the  Original  Values  o7 
Meteorological  Elements 

Strizhevskiy  L.  N.,  Transactions  GMTs,  1968, 

No.  19,  pp.  58-66 

It  Is  assumed  that  observable  values  of  meteorological  elements 
are  the  sum  of  their  actual  values  and  measurement  errors.  The 
statistical  characteristics  of  such  summary  fields  are  known  from 
the  results  of  treatment  of  the  data  of  observations,  and  statistical 
characteristics  of  the  error  fields  were  assigned  on  the  strength  of 
certain  general  considerations  about  the  relationships  of  scales  and 
intensities  of  fields  of  true  values  and  noises. 

On  the  basis  of  the  theory  of  optimum  linear  filtratinn  of  two- 
dimensional  fields  filtration  operators  are  calculated  for  two 
nonlinear  transformations  of  original  fields  widespread  in  meteo¬ 
rology:  the  nonlinear  part  of  vortex  advection  and  the  nonlinear 
part  of  the  individual  wind  derivative.  In  obtaining  the  corre¬ 
sponding  formulas  the  assumption  of  statistical  isotropism  and  homo¬ 
geneity  of  the  original  fields  is  used.  Furthermore,  it  is  assumed 
that  the  investigated  meteorological  elements  are  distributed 
normally..  Calculations  are  conducted  for  various  values  of  scales 
and  intensities  of  the  error  field  corresponding  to  the  various 
degrees  of  complete  treatment  and  the  accuracy  of  measurements  over 
the  different  regions.  The  theoretically  obtained  filtration  oper¬ 
ator  for  the  nonlinear  part  of  the  individual  wind  derivative  is 
used  in  a  one-levei  diagram  of  wind  forecasting,  which  achieved 
definite  Improvements  in  comparison  with  the  variant  of  the  system 
which  did  not  use  smoothing. 

Table  1,  Bibl.  7. 


Objective  Analysis  of  the  Geopotential  Field  Using 
Supplemental  Information 

Sovetova  V.  D. ,  Transactions  GMTs ,  1968, 

No.  19,  pp.  7?-77 

Two  methods  of  objective  analysis  of  a  baric  field  are  described, 
and  results  of  calculations  based  on  one  of  them  are  given. 

The  first  method  is  developed  for  the  geopotential  field  of 
isobarlc  surfaces  of  the  troposphere  (1000,  850,  500  and  300  mbar' 
over  the  northern  hemisphere.  In  poorly  treated  information  of 
regions  of  the  northern  hemisphere  apart  from  data  of  radiosonde 
stations  additional  information  in  the  form  of  results  of  the  mea¬ 
surement  of  pressure  over  the  surface  of  the  earth  was  used:  over 
the  oceans  observations  of  ships  of  the  merchant  marine  were  used. 
Using  mutual  correlation  functions  of  the  geopotential  [7],  the 
method  of  spatial  optimum  Interpolation  [2]  reconstructed  missing 
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information  on  850  mbar  from  the  data  of  measurements  of  the  pressure 
field  over  the  surface  of  the  earth;  at  300  mbar  data  from  500  mbar 
was  used.  Objective  analysis  of  the  geopotential  by  the  method  of 
optimum  interpolation  [3]  was  conducted  at  first  for  one  after  that 
for  another  pair  of  surfaces  with  corresponding  levels  by  the  mean 
autocorrelation  functions.  Results  of  calculations  showed  a  sub¬ 
stantial  improvement  of  the  quality  of  objective  analysis,  especially 
over  oceans . 

The  second  method  was  developed  for  levels  of  the  stratosphere 
over  Europe  and  Asia.  The  problem  of  reconstruction  or  missing 
information  at  high  altitudes  and  subsequent  objective  analysis  of 
the  baric  field  according  to  more  complete  raw  data  of  observations 
was  solved.  As  the  basic  isobaric  surface  according  to  the  data  of 
which  the  information  xas  reconstructed  200  mbar  was  used.  The 
additional  information  in  this  instance  air  temperature  was  used. 
During  the  first  stage  of  solving  the  prob!°">  by  the  method  of 
spatial  optimum  interpolation  the  missing  temperature  measurements 
were  reconstructed  on  100  mbar.  Then  based  on  the  average  tempera¬ 
ture  of  the  layer  and  the  height  of  the  underlying  surface  the  height 
of  the  100  mbar  surface  was  figured.  Objective  analysis  was  carried 
out  always  for  the  same  number  of  stations,  which  made  it  possible 
to  attach  preliminarily  to  each  point  the  eight  best  located  stations 
and  to  interpolate  to  grid  point  according  to  their  information. 


Table  2,  Bibl.  7- 
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IS.  ABSTRACT 

A  finite-difference  model  for  the  solution  of  the  balance  equation 
for  the  northern  hemisphere  is  described.  The  balance  equation, 
rendered  initially  in  a  local  isobaric  system  of  coordinates  (x,  y, 
p)  for  polar  stereographic  projection,  belongs  to  the  family  of  Monge- 
Ampere  equations.  One  of  the  requirements  to  it  is  that  it  must  be 
elliptic  at  all  points  of  the  projection  grid.  The  balance  equation 
was  reduced  to  a  Poisson  equation  with  respect  to  the  function  of 
flux  and  solved  by  the  method  of  successive  approximations.  An 
interesting  characteristic  of  the  solution  is  the  decrease  of  the 
convergence  rate  of  the  iteration  processes  ft  low  latitudes.  The 
field  of  the  flux  function,  obtained  from  the  solution  of  the 
balance  equation  was  used  as  the  initial  field  for  p  quasisolencidal 
forecast,  i,e.,  for  a  barotropic  model.  At  the  end  of  each  forecast, 
i.e,,  every  24,48  and  72  hours,  the  field  of  the  geopotential  is  found 
by  inversion  of  the  balance  equation, i.e.,  by  solving  it  for  the 
geopotential,  using  the  known  distribution  of  the  flux  function.  The 
suggested  method  is  close  to  fast  method  described  by  K.  Mikoyada 
(Numerical  Solution  of  the  Balance  Equation,  Collected  Meteorological 
Papers,  Vol,  X,  No,  1-2,  Tokyo,  i960)  and  by  F.  Schuman  (Numerical 
Methods  in  Weather  Prediction;  The  Balance  Equation,  Monthly 
Weather  Review,  Vol.  85,  No.  10,  1957.  (at  8025208) 
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The  evolution  of  sudden  thermal  perturbations,  which  occur  in  the 
stratosphere  at  the  end  of  winter  and  at  the  beginning  of  the  spring 
and  are  known  as  the  Berlin  Phenomenon,  was  investigated.  It  was 
assumed  that  the  initial  disturbance  of  the  temperature  has  a  spherica: 
shape  and  occupies  a  limited  domain,  and  that  there  is  no  initial 
perturbation  of  the  baric  field.  Finite-difference  equations  were 
approximated,  leading  to  a  numerical  solution,  and  the  solution 
stability  was  tested.  The  resulting  preliminary  data  suggest  the 
following  tentative  conclusions:  1)  The  thermal  perturbations  (at 
maximum  zonal  wind  velocities  of  10  and  30  m/sec)  appear  to  have  no 
substantial  influence  upon  the  atmospheric  circulation,  and 
consequently,  do  not  descroy  the  polar  vortex.  2)  A  perturbation  in 
the  temperature  field  also  causes  a  perturbation  of  the  baric  field 
whose  amplitude  grows  during  the  considered  period  of  time  (4  to  6 
days)  whereas  the  amplitude  of  the  temperature  perturbance  decreases. 
3)  Within  the  first  48  hours,  the  perturbation  maximum  of  the  geo¬ 
potential  remains  at  the  same  latitude  at  which  the  temperature 
perturbation  was  initially  located.  Thereupon,  the  perturbation  of 
the  geopotential  shifts  northward,  whereas  the  temperature  perturba¬ 
tion  remains  at  the  initial  latitude.  This  differs  considerably  from 
processes  that  are  observed  in  nature.  ^1 p  8025209) 
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IS.  ABSTRACT 


Formulas  are  derived  for  the  computation  of  the  weighting  functions 
for  optimal  interpolation  and  smoothing  of  meteorological  elements  for 
the  case  when  the  correlation  functions  of  the  true  (or  discrete) 
magnitude  (signal)  and  of  the  observation  errors  (noise)  are  known. 

For  optimum  interpolation,  it  is  mandatory  to  take  into  account  that 
the  observational  data  contain  random  errors,  and  hence  it  becomes 
necessary  to  obtain  a  magnitude  expressing  the  meteorological  element, 
that  is  free  of  such  errors.  If,  however,  the  errors  are  presumed  to 
include  also  certain  small-scale  fluctuations,  the  assumption  on  the 
uncorrelatability  of  random  errors  would  no  longer  be  valid.  Using 
the  F.  Thompson  technique  for  the  filtration  of  random  processes, 
meteorological  fields  can  be  optimally  smoothed.  The  structure  and 
the  properties  of  the  weighting  functions  are  examined  using  discrete 
as  well  as  continuous  data.  Ultimately,  the  problem  of  determination 
of  a  weighting  function  is  reduced  to  the  solution  of  a  Wiener-Hopf 
integral  equation,  or  to  the  solution  of  a  discrete  analog  of  the  W-H 
equation.  The  W-H  equation  was  solved  using  the  method  of 
indeterminate  coefficients,  approximating  the  correlation  function  by 
a  sum  of  power  functions.  Several  alternate  random  solutions  are 
considered,  e.g.  when  the  errors  are  white  noise,  or  when  the  number 
of  stations  is  increasing  without  restrictions.  Oric  art  has* 

6  figures,  45  formulas.  (AT8025210) 
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An  attempt  was  made  to  derive  formulas  for  the  optimum  smoothing  of 
fields  obtained  from  certain  nonlinear  transformations  of  the  initial 
fields.  The  initial  fields  were  assumed  to  be  isotropic  and  uniform, 
with  a  normal  distribution  and  a  mathematical  expectancy  of  zero.  The 
statistical  characteristics  of  summary  fields  were  derived  from  the 
evaluation  of  empirical  data,  whereas  the  statistical  characteristics 
of  the  error  fields  were  derived  from  certain  general  considerations, 
regarding  the  ratios  of  the  scales  and  intensities  of  the  fields  of 
the  true  values,  and  of  the  noise.  The  filtration  operators  were 
computed  for  two  nonlinear  transforms  of  the  initial  fields,  i.e.,  for 
the  nonlinear  advectlon  component  of  the  vortex,  and  for  the  nonlinear 
component  of  the  individual  derivative  of  the  wind.  The  calculations 
were  performed  for  different  values  of  the  scales  and  intensities  of 
the  field  of  errors  that  corresponded  to  different  degrees  of  exposure 
and  measurement  precision  over  different  areas.  The  theoretically 
obtained  filtration  operator  for  the  nonlinear  component  of  the  wind’s 
derivative  was  used  in  a  one-level  diagram  of  the  wind  forecast.  It 
is  believed  that  this  approach  brings  somewhat  better  results,  as 
compared  to  the  previously  used  variant,  essentially  because  of  the 
more  realistic  values  of  the  wind  velocity  modules.  The  computer 
program  was  also  simplified  and  computer  time  was  reduced.  Grig 
art.  has:  3  figures,  16  formulas,  1  table.  (AT6025211) 
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Two  methods  of  objective  analysis  of  the  baric  field  are  described. 

The  first  method  is  appllcalbe  to  the  field  of  the  geopotential  of 
isobaric  surfaces,  i.e.,  for  four  levels  of  the  troposphere  over  the 
northers  hemisphere.  In  addition  to  the  data  collected  by  radiosonde 
stations  in  areas  of  the  northern  hemisphere,  that  have  a  low  exposure, 
additional  information  was  used:  Pressure  measurements  at  the  surface 
of  the  earth  and  data  collected  by  merchant  marine  vessels.  Using  the 
cross-correlation  functions  of  the  geopotential,  the  information 
lacking  at  the  SpOmb  isobaric  level  was  reconstructed  from  data 
obtained  from  the  measurement  of  the  pressure  field  at  the  earth's 
surface.  The  data  lacking  at  the  300  mb  level  were  reconstructed  from 
measurements  performed  at  the  500  mb  level.  Using  the  method  of 
optimal  interpolation,  the  geopotential  was  analyzed  objectively  first 
for  one  and  then  for  the  other  pair  of  surfaces,  using  the  applicable 
autocorrelation  functions.  The  computations  indicate  a  considerable 
qualitative  improvement  of  the  objective  analysis,  especially  for  the 
water  areas.  The  second  method  was  developed  specifically  for  --he 
stratospheric  levels  of  Europe  and  Asia.  The  problem  involved  the 
reduction  of  information  that  was  not  available  for  higher  levels,  an-: 
a  subsequent  objective  suialysis  of  the  baric  field  from  more  comrlct* 
initial  observation  data.  {AT8#25213} 
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